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Information about the last stages of a binary neutron star inspiral and the final merger can be ex¬ 
tracted from quasi-equilibrium configurations and dynamical evolutions. In this article, we construct 
quasi-equilibrium configurations for different spins, eccentricities, mass ratios, compactnesses, and 
equations of state. For this purpose we employ the SGRID code, which allows us to construct such 
data in previously inaccessible regions of the parameter space. In particular, we consider spinning 
neutron stars in isolation and in binary systems; we incorporate new methods to produce highly 
eccentric and eccentricity reduced data; we present the possibility of computing data for significantly 
unequal-mass binaries with mass ratios g ~ 2; and we create equal-mass binaries with individual 
compactness up to C ~ 0.23. As a proof of principle, we explore the dynamical evolution of three 
new configurations. First, we simulate a q — 2.06 mass ratio which is the highest mass ratio for a 
binary neutron star evolved in numerical relativity to date. We find that mass transfer from the 
companion star sets in a few revolutions before merger and a rest mass of ~ Mq is transferred 
between the two stars. This amount of mass accretion corresponds to ~ 10®^ ergs of accretion 
energy. This configuration also ejects a large amount of material during merger (~ 7.6 x 
imparting a substantial kick to the remnant neutron star. Second, we simulate the first merger 
of a processing binary neutron star. We present the dominant modes of the gravitational waves 
for the processing simulation, where a clear imprint of the precession is visible in the (2,1) mode. 

Finally, we quantify the effect of an eccentricity reduction procedure on the gravitational waveform. 

The procedure improves the waveform quality and should be employed in future precision studies. 
However, one also needs to reduce other errors in the waveforms, notably truncation errors, in order 
for the improvement due to eccentricity reduction to be effective. 
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I. INTRODUCTION 

The majority of coalescing binary neutron star sys¬ 
tems are often expected to have negligible eccentricity, 
low spins, and be very close to equal-mass (with masses 
around l.SSM©). These expectations are based on the 
masses and spins of the population of observed binary 
neutron stars where at least one star is seen as a radio pul¬ 
sar (see, e.g., mm), combined with the efficient shedding 
of eccentricity due to gravitational waves (GWs) during 
the long inspiral mm- However, this observed population 
is quite small, currently consisting of around 12 systems. 
For many of these systems, the evidence that the com¬ 
panion of the pulsar is in fact another neutron star is 
indirect, at best—the companion could still be a fairly 
massive white dwarf (see, e.g.. Sec. 8 in [5] for some dis¬ 
cussion of this issue). Moreover, of these 12 systems, only 
7 have well-determined masses, and only 6 systems (all 
with well-determined masses) will merge within a Hubble 
time, and thus contribute directly to merger rate calcu¬ 
lations; see, e.g., Table 2 of [5] and Table 1 in mm- It 
is thus unclear to what extent these small spins, medium 
masses, and small mass ratios are just a selection effect. 

On the one hand, population synthesis models for bi¬ 
naries formed “m situ” (e.g., 0) predict a much wider 
range of masses and mass ratios than those currently ob¬ 
served; see the discussion in Appendix jA 1| Also the spins 
at merger could potentially be considerably higher than 
those observed in binary pulsars to date, which is dis¬ 
cussed in Appendix [A 2| On the other hand, dynamical 
capture in dense stellar regions, such as globular clusters, 
offers the possibility of forming quite exotic objects, such 
as double millisecond pulsars mm- In fact, there is good 
evidence that dynamical capture and exchange interac¬ 
tions involving neutron stars are a frequent occurrence 
in globular clusters m- Binaries formed by dynami¬ 
cal capture might also have nonnegligible eccentricity at 
merger m- 

Numerical simulations in full general relativistic hydro¬ 
dynamics are the only way to make accurate theoretical 
predictions for the properties of these systems in the time 
period around merger. A prerequisite for all simulations 
are accurate initial data that solve the Einstein constraint 
equations along with the Euler equation on the initial 
hypersurface and also describe the physical system one 
wants to study at some instant of time. Generally, one 
also wants to have this time be not too far from merger, 
to avoid an excessively expensive computation to reach 
the merger. Additionally, quasi-equilibrium sequences of 
initial data at different separations can be used to study 
certain pre-merger properties of these systems without a 
full dynamical evolution. 

Given the potential diversity of the population of coa¬ 
lescing neutron stars in the universe, it is important to be 
able to generate accurate initial data for as much of the 


potential parameter space as possible. In particular, even 
relatively small spin and eccentricity can significantly bias 
measurements of the neutron star tidal deformabilities, 
affecting their ability to constrain the nuclear equation 
of state mils]. Also note that higher mass-ratio neu¬ 
tron star systems, even if rare, are quite interesting from 
a gravitational wave data analysis standpoint, since the 
individual masses of a g = 2.5 system could be measured 
much more precisely than for equal-mass systems m- 

There are a number of well-developed codes for com¬ 
puting binary neutron star initial data in certain portions 
of the parameter space, most notably the open source 
spectral code LORENE [T5]. Other codes include the 
Princeton group’s multigrid solver [16], RAM’s multigrid 
solver |T7], the COGAL code [15], the SpEG code’s spec¬ 
tral solver nMi], and our spectral code SGRID [55H55| . 
All these codes are incapable of reaching certain portions 
of the possible binary neutron star parameter space. In 
particular, they cannot generate consistent initial data 
with (noncorotating) spin with specified eccentricities. 
They also generally have difficulty reaching large com¬ 
pactnesses and high mass ratios. Additionally, when one 
evolves quasicircular initial data computed using these 
codes, one obtains an eccentricity of ~ 10“^, which is or¬ 
ders of magnitude larger than we would expect in most 
binary neutron star systems (see, e.g., i!)- There are 
standard methods for iteratively reducing eccentricity in 
binary black hole initial data (e.g., [55H55] ). However, 
it is not possible to apply these methods to binary neu¬ 
tron star initial data generated with the standard helical 
Killing vector technique. Such methods can only be ap¬ 
plied to consistent initial data if one appropriately gener¬ 
alizes the helical Killing vector, as was only done recently 
in dziiini- Finally, there is no known way of including 
magnetic fields consistently in binary neutron star initial 
data—they are added by hand in all studies that include 
them. 

Of course, it is possible to generate less accurate initial 
data in large regions of parameter space if one does not 
solve the Euler equation and possibly not even the Ein¬ 
stein constraint equations (e.g., using superpositions of 
boosted isolated star solutions, possibly with constraint 
solving, as was done in, e.g., [1511501151] 1. However, if one 
uses inconsistent initial data to initialize a simulation, it 
is always unclear how accurate the resulting simulation 
will be. 

Recently, there has been progress in generating con¬ 
sistent initial data in all these portions of parameter 
space (except for magnetic fields). One of us presented 
a method for constructing consistent binary neutron star 
initial data with arbitrary spin using the constant rota¬ 
tional velocity (CRV) approach in [55] and implemented 
the method in SGRID in [51]. This method has now also 
been implemented by other groups [HHI]. There has 
also been some work on obtaining somewhat high mass 
ratios (up to g = 1.5 [55] [55]) and high compactnesses 
(up to C ~ 0.26 [57]), but neither of these are close to the 
maximum mass ratios (at least ^ 2; possibly up to ~ 3 for 
a large maximum neutron star mass) and compactnesses 
(up to ~ 0.3) that are (at least in principle) possible. We 
also presented a method for generating consistent binary 
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neutron star initial data with arbitrary eccentricity, in¬ 
cluding the possibility of reducing the eccentricity present 
in standard quasicircular data, in m- Concurrently, [29] 
applied a similar method for eccentricity reduction of bi¬ 
nary neutron star initial data. 

Here, we use SGRID to construct binary neutron star 
initial data pushing in all these directions. We have im¬ 
plemented the eccentric neutron star binary initial data 
construction method from m in SGRID, allowing us to 
solve for the velocity potential, which was not feasible in 
our initial multigrid implementation. We have also im¬ 
plemented piecewise polytropes in SGRID, allowing us 
to construct initial data with more realistic equations of 
state (EOSs). We thus show the improvement in the ini¬ 
tial density oscillations of the simple polytropic highly ec¬ 
centric data from m with the new SGRID data, and also 
construct eccentricity-reduced initial data for a simple 
polytropic and a realistic EOS. We also compute aligned 
spin initial data with somewhat larger spins than in |38| 
(as well as for more realistic EOSs). Finally, we illustrate 
the ability of SGRID to compute binary initial data with 
compactnesses up to C ^ 0.23 as well as to compute high 
mass-ratio initial data {q = 2.06) with a realistic EOS. 

Of course, one would like to study the phenomena 
around merger in these newly accessible portions of the 
binary neutron star parameter space: Around merger, 
the strongest gravitational wave, electromagnetic, and 
neutrino emission happens; see |39| for a review of bi¬ 
nary neutron star simulations. Such a study requires 
dynamical evolutions, which will start from initial data 
provided by a member of a quasi-equilibrium sequence. 
In recent years there has been significant work on im¬ 
proving binary neutron star simulations on many fronts, 
notably by including more realistic equations of state 
from piecewise polytropes |40| to finite temperature 
EOSs HU021 , magnetohydrodynamics and neu¬ 

trino cooling HU HU |49] . There are now some simulations 
that include all three of these improvements at once (SO] . 
For our simulations, we use the BAM code in its newest 
version [HHIIKTII^ . 

In this work, we present three different dynamical sim¬ 
ulations. We consider a binary configuration with a mass 
ratio of q = 2.06, which is the largest mass ratio binary 
neutron star system ever evolved in full general relativity. 
This large mass ratio is particularly interesting, since the 
system undergoes mass transfer prior to merger and dur¬ 
ing the merger process a large amount of material gets 
unbound and is ejected from the system. The second 
example we consider is an unequal mass configuration, 
where the two neutron stars have spins misaligned with 
the orbital angular momentum; this is the first simulation 
of the merger of a processing binary neutron star system. 
Here we find the same close relation between the pro¬ 
cessing system (viewed in the nonprecessing frame) and 
the aligned-spin analogue found for binary black holes in 
previous studies (e.g., [53]). As a third test, we perform 
a simulation of an equal mass setup, with and without 
the eccentricity reduction procedure, where one can see 
a clear improvement in the waveform quality due to the 
eccentricity reduction procedure. 

The article is structured as follows: In Sec. ]nj we re¬ 


call the most important equations for our initial data con¬ 
struction and the general framework employed in SGRID. 
In Sec. HI we describe the implementation and the nu¬ 
merical methods focusing on the recent upgrades to the 
code. We summarize the main results of our work in 
Sec. IV where we compute binary neutron star (BNS) 
systems in quasi-equilibrium sequences varying the spin, 
eccentricity, mass ratio, and compactness and show con¬ 
vergence of the numerical method. In Sec.|V| we evolve a 
q = 2.06 nonspinning system, an unequal mass precessing 
configuration, and an equal-mass system with eccentric¬ 
ity reduction. We conclude in Sec. |VI| The appendices 
summarize binary neutron star population synthesis pre¬ 
dictions for more extreme systems, along with some issues 
in predicting the spins expected in binary neutron stars 
at merger; an alternative derivation of the GRV approach; 
and results from our study of single GRV-stars. 

Throughout this work we use geometric units, setting 
c = G = Mq = I, though we will sometimes include 
Mq explicitly or quote values in cgs units for better un¬ 
derstanding. Spatial indices are denoted by Latin letters 
running from 1 to 3 and Greek letters are used for space- 
time indices running from 0 to 3. We always raise and 
lower indices with the physical metric (3-metric for spa¬ 
tial indices and 4-metric for spacetime indices). We shall 
also use index-free notation when convenient, denoting 
vectors (spatial or spacetime) using boldface. 


II. METHOD 


The following section summarizes the fundamental 
framework of the initial data construction in SGRID. 
Since no new development was made in the evolution 
method in this work, we refer the reader to Hi [HU [HIM] 
for details of the methods used in the BAM code. 


A. General Framework 

In this article we investigate BNS systems in quasi¬ 
equilibrium. We construct such configurations with the 
help of a 3 -I-1 decomposition of Einstein’s field equations 
[55] • In this section we recast important equations and 
derive the specific system of partial differential equations 
we solve. 

We start writing the spacetime metric in the form 
ds^ = g^^dx^dx'' 

= —a^dt'^ + 'Yij^dx^ + P^dt){dx^ + /3^dt), (2.1) 

where a is the lapse function and /3* the shift. The spatial 
metric induced on 3-dimensional hypersurfaces of con¬ 
stant t is denoted by "fij. By performing the 3-1-1 decom¬ 
position, the field equations split into two sets, namely 
the Hamiltonian and momentum constraints, 

R- = 167rp, (2.2a) 

Dj = Stt/, (2.2b) 
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and the evolution equations 

— QjOiKij (2.3a) 

dtKij = a{Rij — 2KaK^ i + KKij) — DiDja 

+ CfjK^j — SnaSij + 47rajtj(S — p). (2.3b) 

Here the Ricci tensor Rij and Ricci scalar R are com¬ 
puted from the spatial metric jij with compatible covari¬ 
ant derivative operator Di, and the extrinsic curvature 
Kij is given by 

Kij = ~2^ ~ ’ (2-4) 

where Cp denotes the Lie derivative along the vector field 
f3''. We use K := The source terms (energy 

density, flux, and stress tensor) contained in the right 
hand side of the equations are 


p := 


(2.5a) 

f ■■= 


(2.5b) 



(2.5c) 


where is the normal vector to the hypersurface; we 
also write S := Y^Sij. In this work we assume that the 
matter can be described as a perfect fluid with stress- 
energy tensor 

T>^'^ = [p,{l + e)+p]uV+pg^'', ( 2 . 6 ) 


where po, e, p, and denote the mass density, the in¬ 
ternal energy, the pressure, a nd t he four-velocity of the 
fluid, respectively. Inserting (2.6) in the definition of p, 
j\ and 5*^' [Eqs. (p)5^-([23^ one can obtain explicit 
expressions for the matter quantities entering the right- 
hand sides of the constraint and evolution equations in 
terms of the perfect fluid variables. 

To ensure consistency, we have to solve the constraint 
equations along with the matter equation 


= 0 

and the continuity equation 

S/„{pou'') = 0 , 


(2.7) 


( 2 . 8 ) 


which comes from the conservation of the baryon number. 
Equation (2.7) can be written as the relativistic Euler 
equations 


[po(l -I- e) -bp] u''V^u^ = + u^u'')\7^p. (2.9) 

In many cases it is useful to introduce the specific en¬ 
thalpy 

h = \ -\- e -\- p/ pq. (2.10) 

Then the Euler equations can be written as 

u''V^{hUf,) + Vf,h = 0. ( 2 . 11 ) 

To further simplify the equations, we split the three- 
metric into a conformal factor ip and the corresponding 
conformal metric 7 ij, writing 

Jij = ip‘^%- ( 2 . 12 ) 


Similarly, we express the extrinsic curvature in terms of 
a trace-free piece writing 


K — A ■ ■ 




(2.13) 


Inserting (2.12) and (2.13) in the constraint equations 


gives their final general form before making any addi¬ 
tional assumptions. 


1. Assumptions for metric variables 

In order to obtain a stationary configuration that is ap¬ 
propriate for initial data, we include additional assump¬ 
tions which bring the entire system in an elliptic form. 
The first assumption is the existence of an approximate 
symmetry vector^ 


^ 9 X 9y ') 


(2.14) 


where the functions py are chosen according to the 
problem we want to tackle, allowing us to construct qua¬ 
sicircular, eccentric, and eccentricity-reduced configura¬ 
tions, as discussed in Sec. II B| Together with the exis¬ 
tence of k and the assumption 


Ckg^. = 0, (2.15) 

we assume spatial conformal flatness^ and maximal slic¬ 
ing 


lij = fij, (2.16a) 

K := = 0, (2.16b) 

where fij denotes the flat metric which simplifies to Sij 
in Cartesian coordinates. We preserve these conditions in 
time (at least infinitesimally), that is, dt^ij = Ckjy,i> = 0 
and dtK = CkK = 0. Note that we have not used the 
assumption of an approximate symmetry vector in ob¬ 
taining these last equalities, which is the usual approach, 
with conformal flatness and maximal slicing imposed af¬ 
terwards, as discussed in Sec. Ill A of m- We also obtain 

(2.17) 


^ Notice that in this work we employ a different notation than 
presented in |23II24| where the Killing vector was denoted with 
We change our notation to emphasize that we assume a symmetry 
vector, not necessarily a helical Killing vector. 

^ See, e.g.. Sec. Ill A in El for a discussion of the limitations of 
the conformal flatness assumption, but note that these caveats 
are quite mild for the systems we are considering. In particular, 
while the conformal flatness assumption is a significant obstruc¬ 
tion to constructing high-spin binary black hole initial data (with 
dimensionless spin j > 0.93) as discussed in, e.g., [56] , neutron 
stars cannot spin rapidly enough for this to be a problem. Specifi¬ 
cally, neutron stars have maximum dimensionless spins of at most 
0.7, except for strange quark stars, which can have dimension¬ 
less spins greater than 1; see, e.g.. Figs. 3 and 6 in El- Since 
we are not considering strange quark stars in this work, we do 
not expect that the assumption of conformal flatness places any 
restrictions on the parameter space we can cover. 
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with (L/?)*^ = D'‘— , where Di denotes 

the flat-space covariant derivative. 

Finally, these assumptions lead to the following partial 
differential equations 

= -^{LBYHLBP, - 2TTi^^p (2.18a) 
Dji-LBY^ = {Wy^Dj In -f IdTraV^^j* (2.18b) 

D\ay)=ay (^^yLBr{\.By,+2T:y\p + 2S)^ , 

(2.18c) 

with Di = di in Cartesian coordinates and i3®=/3*-|-fc®-|- 
D.€iji{x^ — x^yYa}, where x^y^ is the center of mass, D. the 
orbital frequency, Ciji is the Levi-Civita symbol, and a* is 
a unit vector pointing along the direction of the orbital 
angular momentum. 


2. Assumptions for matter variables 

Similarly to the metric variables, we also need to make 
assumptions for the matter fields. These assumptions are 
discussed in more detail in [TTl |23] and briefly described 
below. We start by splitting the four-velocity into a piece 
along and one orthogonal to it, which we call V^. 
Specifically, we write 

= u°{k^^ + (2.19) 

with vP = —u^rZjy/a. Next we define 

P/x = hUij,. (2.20) 

While we can assume that CkPfi = 0 for irrotational bina¬ 
ries, this equation is, in general, not satisfied for spinning 
neutron stars—see appendix A of [23) . Thus, we intro¬ 
duce the canonical momentum 1-form of a fluid element 

Pi = h^Pv (2.21) 

and split Pi into an irrotational part which can be written 

as the gradient of a potential. Dip, and a rotational part 

wp. 


Pi = Dip - 1 - Wi 

(2.22) 

or equivalently in four-dimensions 


PiL = Viip + Wfi. 

(2.23) 

Although CkPii 0, we assume 


Ck{pu°) = 0, 

(2.24a) 

l^^Ck{^nP) = 0, 

(2.24b) 

Ai^CkW^ = 0, 

(2.24c) 

with 


IJ 

(2.25) 


which is parallel to the worldline of the star’s center. At 
this point useful relations can be derived immediately. 


W -I- Afc’ 


\2.24h\ 



I'l 

\2.24:c\ 


— 


|2T9) 



AfJ-'kPi^ L,kWu = li' Ck+^k'Utv 

IW,, 




(2.26a) 


_ fc* + Afc* = p ^ 


hu^ 
(2.26b) 


with the three-dimensional Lie derivative and = 
(0, Afc®). Additionally, from the fact that hvP and 'yij are 
approximately constant along we obtain 


(3)£ 


V+Afe 


Wi = 


= +w^ 


hu^ 


r7xi «0. 

(2.27) 

Plu gging ( |2.19|) in to the continuity equation (2.8) and 
using (2.15), ( 2.24a[ ) we get 

A {poau°V^) = 0. (2.28) 


Similar ly, the Eule r equation ( |2.9[ ) together with (2.26aI, 
(2.26b), and (2.27) can be simplified to 


(2.29) 


A [^-,+VW,pj=0, 

which can be integrated to obtain 


—XT -I- Djp = —C = const. 


(2.30) 


Note that a simple derivation of this first integral, which 
makes use of the Cartan identity, can be found in Ap¬ 
pendix 1^ 

The constant C is chosen during the numerical itera¬ 
tion process in such a way that the baryonic mass of each 
star stays constant; see Sec. Ell 
In general the velocity is given by 


which brings the continuity equation in the form 


D,, 


^^{D^P + A) - poau^ip^ + ky 


= 0 . 


(2.31) 


(2.32) 


This equation can be seen as a nonlinear elliptic equation 
in p and especially needs known boundaries at the star’s 
surface to be solved. To handle this issue we introduce 
surface-fitted coordinates in the subsequent section. 
Integrating and using = — 1 leads to 


with 


L^ = 


h = a/A _ _|_ yOi){D^p + wY, 


b + \lb‘^ — 40"^ [{Dip + wpwp 


(2.33) 


(2.34a) 


6 = [(fc*-f/3*)A</> - C] +2a^{DiP + Wi)w\ (2.34b) 


For the data constructed with the CRV-approach we 
choose throughout the entire paper 

A = e*^V(x'=-4j, (2.35) 

where xy^ gives the coordinate position of the center of 
the star and w* is an arbitrary angular velocity vector. 
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B. Specifying the symmetry vector 


Helical Killing vectors are well known and commonly 
used constructs in numerical relativity to construct bina¬ 
ries on circular orbits that are stationary in a corotating 
frame. The general expression for these vectors is given 

by 


= r + 


n, 


qcixy^ - 


where we used the vectors t = dt, x = 


(2.36) 

and 


>'t, X = Ux, y = dy 
if = dip that generate translations in the t, x, and y di¬ 
rections, respectively, and rotations in the tp direction. 
In [17j . we showed how to generalize this vector to in¬ 
corporate eccentricity as well as radial velocity. How¬ 
ever, because our initial numerical implementation of the 
method used a Cartesian grid, without the surface-fitted 
coordinates needed to solve for the velocity potential, we 
settled on a constant fluid-velocity approximation instead 
of solving Eq. (2.321. We are now able to solve the full 
set of equations for the first time. In the following, we 
will briefly summarize how we generalize the standard 
approximate helical Killing vector to an approximate he¬ 
lical symmetry vector that incorporates radial velocity 
and eccentricity. 

To find a vector that approximately Lie-derives the 
flow we make the following two assumptions: (i) Such a 
A:“ exists, (ii) is along the motion of the star center. 

In order to describe eccentric orbits we make the addi¬ 
tional assumption that (iii) each star center moves along 
a segment of an elliptic orbit at apoapsis. Since we only 
need a small segment of an orbit near apoapsis, we will 
approximate this segment by the circle inscribed into the 
elliptical orbit there. Then the radii of the inscribed cir¬ 
cles are 


rc,., = (1 - e)di,2, (2.37) 

where di and d 2 are the distances of the particles from 
the center of mass at apoapsis and e is the eccentricity 
parameter for the elliptic orbit mi- These two inscribed 
circles are not centered on the center of mass, but on the 
points 

= ^ 1,2 T ?’ci,2 = a^CM + e(a;i,2 - sjcm), (2.38) 

where we have used di ^2 = |a:i .2 — a;cM| and assumed 
that apoapis occurs on the a;-axis. (The upper and lower 
signs correspond to the subscripts 1 and 2, respectively.) 
Assumption (ii) then tells us that the approximate Killing 
vector for elliptic orbits must have the form 

fc“cci .2 =t°‘ + ^[{x- a;ci ,2 )y“ - y a;“] (2.39) 

near each star m- 

The next step is to allow a slow inspiral of the orbit due 
to energy loss because of GW emission. This means the 
orbital velocity will have a small radial component in the 
direction of the center of mass. Assumption (ii) then tells 
us that assumption (iii) above needs to be modified to 
include a radial piece. We assume that the approximate 
Killing vector now is 

ki,2 = Kcci,2+—r°‘ = t“-fH[(a;-a;ci,j2/“-?/a;“]-f—r“, 

ri2 ri2 

(2.40) 


which we also refer to as a helliptical approximate sym¬ 
metry vector. Here r“ = (0,x,y,z) points in the radial 
direction, ri 2 = |a;i —X 2 \ is the distance between the star 
centers, and Vr is a radial velocity parameter. The radial 
velocity Vr could be chosen corresponding to the radial 
velocity of an inspiralling binary from post-Newtonian 
calculations, or it can be obtained from an iterative pro¬ 
cedure aimed at reducing the orbital eccentricity such as 
the one described in Sec. lIV Cl In this case we also have to 
adjust the eccentricity parameter e that appears in Xci 2 - 
The reason is that changing e amounts to changing the 
tangential orbital velocity, which is needed when we want 
non-eccentric inspiral orbits. 

To have a more physical quantity that will be useful 
for comparisons, we consider the mean motion, as in E!: 

n = 2 tt/T = H(1 -k e) Vl - e2, (2.41) 

where T is the orbital period. 


III. CODE DESCRIPTION 


To construct initial data with SGRID we use the nu¬ 
merical framework described in El 121] • In this section 
we recall important aspects to give an almost complete 
picture. In particular, we describe the grid configuration, 
the iteration procedure, and recent changes which allow 
us to use more realistic EOS and also compute configura¬ 
tions with a relatively large mass ratio. A short discussion 
about the BAM code will be given in Sec. 


A. Grid configuration 


We place the neutron stars along the x-axis. Eigure[^ 
illustrates the part of the computational domain with y > 
0, z = 0. The entire grid is built out of six individual 
domains. The grid configuration is not fixed and changes 
during the computation in response to changes in the 
positions of the stars’ surfaces. 

We follow the approach of [22l [58] and express Garte- 
sian coordinates as 

^ 2 ((A2 + R2)2 + l) (3-la) 

y = ^( (A2+R2)2 (3-lb) 

^ = ^( (A 2 + R 2)2 +i) XRsiniip), (3.1c) 

with X G [0,1], R G [0, -s/l — A2], tp G [0, 27r). Further¬ 
more, we transform to coordinates A, B, ip for the differ¬ 
ent domains. The domains covering the exterior of the 
stars and including spatial infinity (A, B) = (1, 0) employ 
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A=1 


X 


A=1 


FIG. 1: The grid structure in the a:t/-plane for an equal-mass configuration. One can see lines of constant A and B (dark blue 
lines) for b = 16, (t+ = —a- = 1.304 and jdmax = 0.5 with an overlay of the density profile of the cross section. Moreover one 
can see the Cartesian boxes with Chebychev grids inside the stars. 


with 


X = il-A) [Re(C'± {B, ^)) - B Re(C± (1, <^))] 

+ Bcos + (1 - ^)arg(C'±(l,y)))^ , {3.2a) 

R={1- A) [lm{C±{B,ip)) - BIm(C±(!,</?))] 

+ Bsin + (1 - A)arg(C±(l,<p))^ . (3.2b) 


D± := {I - A) aig{C±{l,ip)) + 5^'^A, (3.7) 

where (5± = 1 for the star with a; > 0 and is zero for the 
other, and vice versa for i5;p. Unfortunately, the trans¬ 
formation to (^, B, ip) coordinates is singular for A = 1 
(i.e., at the star’s center). To cure this problem, we cover 
the center by a Cartesian box with grid points at 


Since spatial infinity is included in our domain, we can 
impose exact Dirichlet boundary conditions: 


X k = 


-cos 


kn 
d - 1 


+ 


(3.8) 


lim ij) 

= 1, 

(3.3a) 

lim B* 

r—>-oo 

= 0, 

(3.3b) 

lim aij} 

= 1, 

(3.3c) 


r—^oo 


where r denotes the coordinate distance from the origin. 

The inner domain boundary A = 0 is the star’s surface, 
with 


C±{B, if) = Wtanh 


^±{B^ + t7rB\ 


■ 


(3,4) 


where (j± is a function that determines the shape of the 
star’s surface, and ± denotes the sign of the ^-coordinate, 
i.e., the left or the right star. At each star’s surface 
Eq. (2.32[) is subject to the boundary conditions 


[(B* -b wV) - -b fc*)] D,p = 0. (3.5) 


where x* = {x,y,z), with 0 < k < ncart- The Cartesian 
boxes cover a region for A > Amax- The choice of A^ax 
allows us to specify the clustering of the grid points. For 
large A^ax the Cartesian box is smaller, while for small 
Amax the box is larger. Thus, introducing a small A^ax 
increases the resolution in the outer region of the stars. 
This will be important when piecewise polytropes are em¬ 
ployed, where it is crucial to resolve the crust with a suf¬ 
ficient number of grid points. The collocation points in 
the other regions of the grid are 


Ai 

Vk 


A„ 


2 

1 
2 
2'Kk 


1 — cos 


TTl 


riA - 1 


1 — cos 


TTJ 


ns - 1 


(3.9a) 

(3.9b) 

(3.9c) 


The coordinate transformations inside the stars are 

A = (1 - A) [Re(C±(B, p)) - B Re(C±(1, (^))] 

-bBcos(B±)-b(5±(l-B)A, (3.6a) 

B = (1 - A) [Im(C'± (B, if)) - B Im(C± (1, tp))] 

-bBsin(B±)-b5=F(l-B)A, (3.6b) 


with 0 < i < HA, 0 ^ j < nB, 0 < (f < n^. In the 
A, B-directions we use Chebyshev polynomials, while for 
the (/J-direction a Fourier expansion is used. For a typical 
configuration, we employ between 20 and 28 points in 
A, B and 8 points in the (/3-direction. The Cartesian box is 
covered with = Uy = = ncart = 16,..., 24 points, 

where typically we choose ncart = — 4. 





































equation (3.16); see IIIC2 for more details. 


(vi) We then compute h and choose C± such that 
the baryonic mass of each star remains constant. 
Afterwards, we update a± to reflect the changes in 
the shape of the stars’ surfaces and adjust the do¬ 
main boundaries accordingly. In most cases we filter 
out high frequencies in a± for overall stability and ap¬ 
ply dBO'±{B, ^)\b=o,i = 0 to keep the stars on the x-axis. 


(vii) We go back to step (ii). 


C. Code Improvements 


FIG. 2: Iteration scheme as outlined in the text. 


1. Including piecewise polytropes 


Finally, some regularity conditions along the x-axis 
have to be imposed: In the domains where (A, B, cp)- 
coordinates are employed, we set 

d^F = 0, (3.10a) 

dsF + dsd^d^F = 0, (3.10b) 

with F G {tp,B\a,4>} and s := 


B. Iteration procedure 


To solve the coupled system of partial differential 
equations we perform a specific iteration procedure first 
introduced and described in detail in [24]. The scheme 
is sketched in Fig. and we describe it in detail in the 
following: 


(i) We start with an initial guess. This guess can be 
obtained by solving the Tolman-Oppenheimer-Volkoff 
(TOV) jSni ISD] equation (and superposing, if we are con¬ 
sidering a binary) or given by a previously constructed 
configuration. The velocity potential (j) in each star is set 
to 0 = Q.{xc* — xcM)y-i where xc* is the cc-coordinate 
of the star’s center (this initial guess corresponds to a 
spatially constant velocity field and is exact for rigidly 
rotating nonrelativistic stars). 


(ii) In the second step we evaluate the residuals of all 
elliptic equations (denoted by Aeu in Fig. and stop if 
these residuals are below the prescribed tolerance. 


(iii) If the residual of Eq. (2.32) is bigger than the 
combined residuals of Eqs. (2.18), we solve (2.32) for cj) 
and use a softening procedure (p = C</'soived + (I — C)?^oidj 
where for this iteration procedure ( = 0.2 is applied. 


(iv) We solve the elliptic equations for '0,1?*, a (2.18) 
with a softening of C = 0.4. 


(v) The positions of the stars’ centers, xc*,±, are 
determined by the maximum of h along the cc-axis. We 
determine and xcm with the help of the force balance 


In order to easily incorporate more realistic equations 
of state we follow the approach in m and approximate 
them by piecewise polytropic equations of state. 

For a simple polytropic equation of state (EOS), with 
p = kpq (where F = I -|- I/n), the matter-variables inside 
the star are C°° and only the star’s surface needs spe¬ 
cial attention. However, this is no longer the case when 
dealing with piecewise polytropes. 

For a polytrope with polytropic index nj and constant 
Kj, the pressure p, specific enthalpy h, and energy density 
Pe are related to the rest mass or baryonic mass density 
Po via 


p = 

l + l/ni 

KIPO 

(3.IIa) 

h = 

{ni + 1 )kipI^'^^ + Kj, 

(3.IIb) 

Pe = 

{niKipp"“ + Ki)po, 

(3.IIc) 


where Kj is a constant that determines the specific en¬ 
thalpy at the star surface. For piecewise polytropic EOSs 
we divide the range of possible po into intervals [0,po,i], 
[Po,i)Po. 2 ]) etc. We label the intervals by J = 0,1,.... 
Within each interval I we use the polytropic relations 
of Eq. (3.1 la), but with a different nj, kj and Kj. In 
interval 0 we must choose Kq = I to ensure that the spe¬ 
cific enthalpy is unity at the star surface. We can freely 
choose all the nj and kq to closely approximate some de¬ 
sired EOS. However, in order to assure continuity of p, h 
and Pe across interval boundaries the remaining ki and 
Ki must be related by 


Kl 

Ko 

Ki 




1 , 

Ki^i + m-iKi^ipl^j 


(3.12a) 

(3.12b) 

njKjpl'^plS.Uc) 


In this case we only know that all the matter variables 
are at least C° (i.e., continuous), but not necessarily dif¬ 
ferentiable.^ We use the parameters presented in m 


^ We have not used the spectral fits to realistic equations of state 
from m, even though these would give matter variables 
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FIG. 3: The mass-radius relation for the EOSs employed in 
this work. See Table U for further details. 


and employ four different pieces consisting of a crust and 
three inner regions. All piecewise polytropic EOSs we use 
give a maximum mass of Mmax > I.QQMq (so they are 
la compatible with the precise high-mass neutron star 
measurements in |631164j l. and have an adiabatic sound 
speed Cg < 1 for densities up to the maximum density 
of a stable TOV star. These EOSs also span a range 
of microphysical content, including hyperons (H4) and 
the hadron-quark mixed phase (ALE2), as well as the 
standard npe/r composition, and were obtained using a 
variety of calculational methods (see Sec. II in [ST] for 
further details). Important parameters for the EOSs we 
employ are given in Tab.jTjand the mass-radius relations 
for TOV stars are shown in Fig. [^ Additionally, we use 
two simple polytropic EOS, one (r2) with n = 123.6489 
and r = 2 and one (r2.72) with k = 23841.43 and 
r = 2.7203. The latter EOS is constructed as an “av¬ 
erage fit” to certain of the realistic EOSs we consider, 
where we have fitted p{po) = PiiPo)/^ for the six 

EOS S = (SLy, ALF2, MPAl, MSI, MSlb, ENG) with 
a simple polytropic EOS in the density interval po G 
[0,3.24 X 10“^]. The resulting EOS tends to be more 
“realistic” than the r2 polytrope widely used in the liter¬ 
ature. In particular, it allows a maximum mass > 2Mq 
and has a maximum adiabatic sound speed < 1 for the 
maximum mass TOV star. The overall qualitative behav¬ 
ior of its mass-radius curve is also more similar to those 
of realistic EOSs than that of the r2 EOS (see Fig. |^. 

In the past, SGRID used q := p/po as the fundamen¬ 
tal variable, i.e., the matter variables and their spatial 
derivatives were all derived from q. But with this defini¬ 
tion q will only be in case of a piecewise polytropic 
EOS, even for a single TOV star. In contrast, h, which 


TABLE I: Properties of the equations of state (EOSs) used in 
this work. The first seven rows refer to piecewise polytropes, 
where we employ the fits of m- These EOS use a crust with 
Kcrust = Ko = 8.948185x10“^ and Ecrust = l-|-l/no = 1.35692. 
The divisions for the individual parts are at po^i = pcrust, 
po ,2 = 8.12123 X 10“^ and po ,3 = 1.62040 x 10“^. The last 
two rows refer to simple polytropic EOSs. The columns (for 
the piecewise polytropes) refer to: the name of the EOS, the 
maximum density in the crust, the three polytropic expo¬ 
nents T/ = l + l/nr for the individual pieces, and the maxi¬ 
mum supported gravitational mass maximum baryonic 

mass and maximum compactness respectively, of 

an isolated nonrotating star. [We define the compactness by 
C ~ M/R, where R is the star’s radius (in Schwarzschild 
coordinates) and M is its gravitational mass.] For the sim¬ 
ple polytropes we present F and k, in addition to the same 
maximum values for an isolated nonrotating star given for the 
piecewise polytropes. 


EOS 

pcrust ' 10 

Fi 

F2 

Fa 


jy^max 

^max 

SLy 

2.36953 

3.005 

2.988 

2.851 

2.06 

2.46 

0.31 

ALF2 

3.15606 

4.070 

2.411 

1.890 

1.99 

2.32 

0.26 

ENG 

2.99450 

3.514 

3.130 

3.168 

2.25 

2.73 

0.32 

H4 

1.43830 

2.909 

2.246 

2.144 

2.03 

2.33 

0.26 

MPAl 

2.71930 

3.446 

3.572 

2.887 

2.47 

3.04 

0.32 

MSI 

1.52594 

3.224 

3.033 

1.325 

2.77 

3.35 

0.31 

MSlb 

1.84169 

3.456 

3.011 

1.425 

2.76 

3.35 

0.31 

r2 

r = 2 . 

n = 

123.6489 

1.82 

2.00 

0.21 

r2.72 

F = 2.7203, K = 

23841.43 

2.40 

2.85 

0.30 


is given by 


dh ^ h[m{r) + 47rr^p] , 

dr r[r —2m(r)] 

for a single TOV star, will be at least inside the star 
under the assumption that p,po G C®. [Here m(r) := 
47r po(f)f^df, so m € .] For this reason we have 
switched variables to 


q:=h-l, (3.14) 


which is at least for single TOV stars. Taking spatial 
derivatives of q is thus more accurate than taking them 
of p/pq. 

We can compute the other matter variables in terms of 
q, giving 


Po 

P 

Pe 


+ l-Ki 

_Ki{ni + 1 )_ 

q + l-Ki 

Po -—j 

ni + I 

nip + KiPq. 


(3.15a) 

(3.15b) 

(3.15c) 


inside the star, since we need to have the same implementation 
of the EOS as in the BAM code that we use for evolutions: Slight 
differences between the two fits to a given EOS would lead to 
unphysical effects upon starting the evolution. 


2. Updating and 


In order to also solve Eq. (2.331 we need to know the 


values of Vt and We first determine the star centers 
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^C*± finding the maximum of the current h along the 


x-axis. Using dih\^i^ ^ = 0 in Eq. (2.33) we find 


di In 


w 




(a 


, '^2 \ 

+ 


= -29ilnr| , (3.16) 


Note that + A:* is a function of fl and The right 
hand side of Eq. (3.16) is given by 


au 


r = 


1 - 


(/3* 




Djcj) 


hu^ j OL^hu^ (a/iu°)2 


^1 - (^/ 3 * + fc* + (Pi + ki + ^ 


, (3.17) 


which we reproduce here (with replaced by fc*) because 
there were some typos in the original published version 
in [H]. Eq. (3.16) is called force balance equation. It 
gives one equation for fl and at each star center and 
thus can be used to update fl and One notewor¬ 

thy cavea t is that we evaluate the derivative of InE in 
Eq. (3.16) for the O and before the update. 

We have found that the force balance equation works 
well in many cases. However, for more massive stars or 
higher mass ratios the overall iteration can become un¬ 
stable. In this case the center of mass drifts away and the 
magnitude of the Arnowitt-Deser-Misner (ADM) momen¬ 
tum 


P: 


ADM 


= / iV 




X, 


(3.18) 


especially of its j/-component T^dm becomes large. This 
problem has also been observed by others m- It can 
be solved in part by computing fl and in a different 
way: Notice first that the matter flux 

f = a{pE + p){u°)\V^ + e+ /?*) (3.19) 


in Eq. ( 3.18[ ) depends on the Killing vector and thus 
on H and~EQ]y[. Using the D from the previous iteration 
we can then solve the equation P^^m ~ b for This 

gives a value for the center of mass such that the P^um 
will be zero as desired. Once we have determined a::cM 
in this way, we next compute U from Eq. (3.16) for each 


star’s center. This will in general give two different values 
for U. For the final SI we simply use the average of these 
two values. 


IV. BINARY NEUTRON STARS IN 
QUASI-EQUILIBRIUM 



FIG. 4: Reduced binding energy vs. specific total angular mo¬ 
mentum for a binary system with = Mjf = 1.4895 and 
the SLy EOS. The influence of aligned and antialigned spin is 
shown with arrows 



e 


FIG. 5: The spin-orbit coefficient £so in the reduced bind¬ 
ing energy of the binary neutron stars. We compare numer¬ 
ical results for the SLy (top panel), r2 (middle panel), and 
MSlb (bottom panel) EOSs with the predictions of the post- 
Newtonian (PN) approximation (black dashed lines). We also 
include an average line for our data in all panels (solid black 
lines). We estimate the numerical error bars from computa¬ 
tions at different resolutions and show them as shaded regions. 


In the following, we discuss the main results on equilib¬ 
rium configurations obtained by applying the framework 


and code improvements discussed in Sec. [H] and Sec. Ill 


In particular we analyze the spin-orbit (SO) interaction 
for realistic EOSs; we extend the work of m on highly 
eccentric orbits, where we improve our data by solving 
Eq. (2.32) for the velocity potential; we investigate inspi¬ 
rals on eccentricity reduced orbits; we also present signifi¬ 
cantly unequal mass setups as well as configurations with 
high compactnesses; we end with a convergence study. 


A. Spins 

To date, there exist three ways to model the fluid ve¬ 
locity field in quasi-equilibrium BNS configurations in 
general relativity. Neutron star spins are neglected in 
the irrotational approach, while they are treated in an 
unphysical manner if corotation is assumed. The CRV 
approach [521 [21] (see also [25|) is the only known alter¬ 
native for the construction of consistent and constraint 
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solved initial data where spins are included and can be 
chosen freely. Alternative approaches have been proposed 
and employed in dynamical evolutions [311138 j . but they 
are based on constraint-violating data which also violate 
hydrodynamical stationarity conditions. Here we discuss 
some properties of equilibrium sequences of CRV BNS ex¬ 
panding the work of [38], and focusing on the spin-orbit 
interaction. 

We compute equilibrium sequences for the SLy and 
MS lb EOSs (Table [I]), setting the baryonic masses of the 
individual neutron stars to Mi, = 1.48945. Additionally, 
we use results of [38], where the T2 EOS was employed 
and the baryonic mass of the neutron stars was set to 
Ml, = 1.625. For each of these EOSs we obtain sequences 
at fixed baryonic mass and for five different spin magni¬ 
tudes, two aligned spin setups (ft)) one irrotational (00), 
and two anti-aligned setups (ii). Each sequence essen¬ 
tially mimics an adiabatic evolution. (In the aligned and 
anti-aligned cases we are considering, the spin directions 
remain unchanged during an evolution.) The resolution 
employed for the piecewise polytropes is = 28; 

riy, = 8, ncart = 24, while for the r2 simple polytrope it is 
riA = riB = M, = 8, ncart = 20. We use a higher res¬ 
olution for the piecewise polytropic r uns to better resolve 
the crust region (cf. Table [I] and Sec. IV F). 

We stress that an unambiguous definition of the indi¬ 
vidual spins of the stars in a binary system is in general 
not possible. It is however possible to define the spin of 
single isolated neutron stars within the CRV approach, 
as discussed in Appendix 0 Using this result, one can 
give an estimate of the magnitudes of the spins of the 
stars in a binary by considering, for each star, an isolated 
configuration with the same baryonic mass Mi, and the 
same rotational part of the 4-velocity w^, and assuming 
that the individual spin magnitudes S^, remain un¬ 
changed when we use the same parameters to compute 
binary initial data. 

We analyze equilibrium sequences in terms of the re¬ 
duced binding energy 


^ _ 1 / Madm ^ 


(4.1) 


and the specific orbital angular momentum 


i = 




Jadm -S^-S^ 
uM^ 


(4.2) 


Here i/ := M^M^/M"^ is the symmetric mass ratio, M^ 
and M^ are the individual masses of the stars in iso¬ 
lation, M := M^ + M^, and Madm and Jadm are 
the binary’s ADM mass and angular momentum, respec¬ 
tively. Here E(i) is a gauge-invariant way to characterize 
the dynamics, which is also applicable to full numerical 
relativity evolutions; see [SS]]^ and Sec. VB The E{i) 
curves are shown in Fig. [^for the SLy EOS; other EOSs 


Notice that the individual masses of the stars in isolation (M^, 
M®) are obtained here for spinning neutron stars and differ from 
the results for irrotational stars with the same baryonic mass. 


show qualitatively the same behavior. From the hgure 
one observes that aligned configurations are less bound 
than antialigned configurations. This behavior follows 
from the fact that the spin-orbit (SO) interaction, which 
is the main spin-related effect, is repulsive for aligned 
spins and attractive for antialigned spins (for a hxed £), 
see, e.g., [67] and below. 

In the following we explicitly compute the SO contri¬ 
bution in our binding energy data, and show its influence 
on the dynamics as well as its independence of finite size 
(EOS) effects. We write the binding energy as 

E{£) = Eq + Eso + Et + Ess , (4.3) 

where Eq describes the binding energy of a nonspinning 
black hole binary in the conformal flatness approxima¬ 
tion, and is therefore independent of the spin and mat¬ 
ter effects; Eso and Ess represent the SO and spin-spin 
contributions, respectively; and Et denotes the tidal con¬ 
tribution. Assuming for simplicity that the dimension¬ 
less spins are the same (as is the case here) j := = 

/{M^Y = 3^ = /{M^Y, the SO interaction is 

proportional to 

Eso^ 3 -L = \\j\\\\L\\ cos(Z(j, Lj). (4.4) 

Thus, the angle between j and the orbital angular mo¬ 
mentum L dehnes whether the SO interaction is repulsive 
or attractive. In the cases considered here, cos(Z(j,L)) 
takes the values 1 and —1. For these two possibilities we 
can write 


Eso=£soie)j, (4.5) 


where j denotes the signe d m agnitude of j, i.e., j := 
IIjII cos(Z(j,L)). Equation ( [4.5] ) allows us to answer two 
important questions: (i) Do we see an imprint of the EOS 
on the SO interaction? (ii) Does the linear dependence 
of Eso on j capture the main dynamics? 

According to our spin definition (see above), the spins 
are constant during the adiabatic evolution. This is a 
good approximation at these separations and also sup¬ 
ported by numerical evidence in binary black hole simu¬ 
lations [55H7T] . The £so{£) term can be computed using 


£so{^) = 


Ef^\£)-Ef^\£) 

*2j 


(4.6) 


for different (signed) spin magnitudes j, where E^^^\£) 
and E^^^\£) is E{t) for aligned and anti-aligned spins. 


Indeed, for a given spin, all the terms in (4.3) except Eso 


cancel in the combination (4.6) because they all have the 
same sign. 

The function £so{^) is shown in Fig. for all three 
EOSs considered here: SLy (top panel)) r2 (middle 
panel), and MS lb (bottom panel). Additionally, we com¬ 
pute the average for all EOSs (solid line) and compare 
our results with the linear-in-spin part of the 4PN en¬ 
ergy from Eq. (8.23) in [75], which is shown as a black 
dashed line in Fig. From the figure we observe: (i) 
£so{£) is positive, therefore the SO-interaction is repul¬ 
sive/attractive (positive/negative) according to the sign 
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FIG. 6: Star’s maximum density scaled by its initial value: 
Simple superposition of TOV stars without solving the con¬ 
straints (dashed green), Cartesian multigrid method solver 
of |17| (solid red), SGRID data assuming constant three- 
velocity (dotted blue) and SGRID data solving the fluid po¬ 
tential equation (|2.32[) (solid black). 


of j [in general, this depends on Z(j, L)]; (ii) all the curves 
agree within their errors, i.e., there is no significant de¬ 
pendence on finite size (EOS) effects, (iii) the PN expres¬ 
sion from captures the behavior of our conformally 
flat data for all employed EOSs. 

We notice that the Ess and the Ej- terms in (4.3) 
can be extracted in a similar way. However, for the spin 
magnitudes and orbital separations considered here they 
lie within the uncertainty of our data. 


B. Highly eccentric configurations 

In [T7] we described a method to produce hydrody- 
namically consistent initial data for relativistic stars on 
orbits with arbitrary eccentricities for the first time. In 
our initial implementation of the method, we used an el¬ 
liptic solver based on a Cartesian multigrid method, for 
which it is technically difficult to solve an equation to 
determine the velocity potential, since this would require 
boundaries at the star’s surface. Instead of introducing 
surface-fitted coordinates, we employed a constant three- 
velocity approximation, which could be motivated by the 
restriction to irrotational binaries. This means that we 
assumed the instantaneous (at apoastron) three-velocity 
of a fluid element measured by a coordinate observer 
to be constant throughout the star, so the four-velocity 
could be written as 

M“ = M*(r+z;^j/“). (4.7) 

However, SGRID provides surface-fitted coordinates and 
allows us to solve easily for the velocity potential. Eirst, 
we want to briefly compare our old results to the newly 
obtained SGRID results and show the improvement of 
the initial data gained by solving the additional equation 
for the velocity potential. Fig. compares the oscilla¬ 
tions of the central density throughout the evolution for 
the previous multigrid solver and SGRID for two equal 


mass stars with baryonic masses = 1.620 on a qua¬ 

sicircular orbit with an initial (2, 2) mode GW-frequency 
of Ma ;^2 = 0.053. 

This simple test case with two polytropic stars (r2) 
clearly shows the influence of the fully solved velocity po¬ 
tential: While we observe strong oscillations of 30% in the 
central density for superimposed TOV stars (which gives 
constraint violating initial data), we only observe roughly 
4% oscillations for the constraint solved data using the 
constant 3-velocity approximation. In this case, SGRID 
and the Gartesian multigrid data give a good agreement. 
If we drop the approximation and solve for (j), we can ob¬ 
tain even lower oscillations, improved by a factor of five, 
i.e., less than 1% (solid black line in Fig. [^. Here we 
use the same evolution setup in BAM as in ^7j , in order 
to make a direct comparison: We use the Baumgarte- 
Shapiro-Shibata-Nakamura IMS] formulation, and 98 
points in each direction in each of 5 refinement levels, 
with a grid spacing of = 0.1875 on the finest level. We 
do not use the conservative mesh refinement introduced 
in [55] ■ The SGRID initial data use = ub = 24, 
= 8 and ncart = 20 points. The reduced density oscil¬ 
lation with our new setup will allow us to study orbit in¬ 
duced oscillations as in [30] in more detail and disentangle 
the orbital effect from oscillations due to the initial data 
construction, which were present in previous attempts. 

We now consider quasi-equilibrium sequences for differ¬ 
ent eccentricities. We compare our eccentric sequences to 
post-Newtonian results for the eccentricity in Fig. [7] Here 
we compute the PN eccentricity from the ADM expres¬ 
sions for the energy and angular momentum; we subtract 
the gravitational masses of the stars in isolation from the 
total ADM energy to obtain the binding energy that en¬ 
ters the PN calculation. (We did not perform such a 
comparison in [13, since the BAM implementation does 
not use a compactifled grid including spatial infinity and 
was thus unable to compute the ADM quantities suffi¬ 
ciently accurately for this comparison.) In this figure, we 
see that the PN eccentricity indeed converges nicely to 
the input eccentricity as the binary’s separation increases, 
but it only converges well to the input eccentricity as one 
increases the PN order for smaller eccentricities (0 and 
0.1). For the two larger eccentricities we consider (0.4 
and 0.5), the IPN results are closer to the input eccen¬ 
tricity than the 3PN results (though the 3PN results are 
closer than the Newtonian or 2PN results). 

In order to obtain the PN expressions we used to cre¬ 
ate this figure, we, in essence, derived the 3PN exten¬ 
sion of the IPN expression for the eccentricity given in 
Eq. (2.36) of Mora and Will [76], using the general 3PN 
results they give. (The only difference in the derivations 
is that since we solve for the square of the eccentricity, 
we do not perform the PN expansion of the square root 
given in their expression.) We start from the expressions 
for E := E},/{Mv) and J := (the scaled binding 

energy and angular momentum) that Mora and Will give 
in Eqs. (2.35) and (2.40) for both harmonic and ADM 
coordinates. (Note that E and J are the same as the E 
and I used in Sec. [IV A[ the latter agreement only holds 
for the zero spin case we consider here.) We then invert 
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FIG. 7: The two plots in the top panel show the squared eccentricity calculated from the ADM energy E and angular 

momentum J using harmonic coordinate PN expressions. We show the results at 3PN (top left) and IPN (top right) for runs 
with different input eccentricity e (given as horizontal dashed lines). (We plot the square of the PN eccentricity since this 
quantity becomes negative for small e.) While the results for small eccentricities are calculated more accurately with increasing 
PN orders, the error is no longer monotonic with PN order for larger eccentricities. The two plots on the bottom panel show 
convergence of the eccentricity estimate with increasing PN order for an input eccentricity e = 0.1 on the left and how this 
convergence becomes more erratic for higher eccentricities (using e = 0.5 as an example) on the right side. 


the series for E to express the PN parameter C, in terms 
of E] we now use E as our expansion parameter. We 
can thus calculate the series for EJ^ by substituting the 
series for ( in terms of E into the PN series for J and 
expanding consistently. 

We then use the resulting harmonic coordinate expres¬ 
sion to obtain a PN series for the square of the eccentric¬ 


ity parameter, in terms of E and J. We exper¬ 

imented with various treatments of the expansion (e.g., 
not expanding after substituting the expression for ^ in 
terms of E and/or solving for the eccentricity numerically 
instead of by series inversion) and found that a consis¬ 
tent expansion to a given PN order produced the most 
reasonable-looking results. Specifically, this gives 


= l-2C+[-A-2rj+ (-1 -b 3r])^]E + 


20 - 23r] 


- 22 -b 60r) + irf - (Sir? -b 


E^ 


-b 


-2016 -b (5644 - 1237r2)7? - 2b2ri^ 4848 -b (-21128 -b 3697r2)77 -b 298877^ 


12^2 


-b ( -30?7-b + ^ 


E-^ 


24 ^ 


- 20 -b 29877 - I 86772 - 477 '^ 


(4.8) 


where ^ := —EJ^. We also performed the same calcula¬ 
tion with the ADM coordinate expressions, which differ 
from the harmonic coordinate ones starting at 2PN, and 
found the expected small differences in the results (< 5%, 
with smaller differences for larger eccentricities). 

Our original motivation for using these particular PN 
results was to avoid using expressions which are not well- 


behaved in the limit of a head-on collision from rest, 
which is what our eccentric data approach in the limit 
e ^ 1, e.g., expressions that have factors of J in the 
denominator. However, even though we start from ex¬ 
pressions that are well-behaved in this limit (by consid¬ 
ering EJ^ instead of J by itself), we still obtain fac¬ 
tors of 1/^ in the final expression for [Eq. (4.8|] 
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starting at 2PN. This is not a significant problem since 
•/adm is not very small for the systems we are consid¬ 
ering. We thus also investigated the expressions given 
by Memmesheimer, Gopakumar, and Schafer (MGS) [77] 
which are not well-behaved in the limit of constant energy 
and vanishing angular momentum. Here one can com¬ 
pute the three eccentricities defined in the post-Keplerian 
parametrization of an eccentric orbit (e^., e^, and e^) from 
Eqs. (20) and (25) in MGS (these give the expressions in 
ADM and harmonic coordinates, respectively). 

We find that (particularly in harmonic coordinates) 
agrees quite well with the input eccentricities, especially 
at large separations, with fractional errors of < 2.4% 
for all separations we consider for e = 0.5. In addition, 
we can compute the coordinate separation r of the stars 
from the ADM energy and angular momentum by not¬ 
ing that the binary is at apoapsis, so we can take r = 0 
in Eqs. (Al) and (A3) in MGS (again, these give the 
expressions in ADM and harmonic coordinates, respec¬ 
tively) and then solve for r numerically. Here we find 
that the value for r we obtain by solving the harmonic 
coordinate equation agrees quite closely with the coordi¬ 
nate separation from SGRID for all the separations we 
consider, particularly for the two larger eccentricities we 
consider (fractional errors of < 2%). Such close agree¬ 
ment may simply be fortuitous, since there is no reason 
a priori to expect the coordinate systems used to agree 
so closely. 


C. Orbits with reduced eccentricity 


The helliptical approximate symmetry vector (or, more 
formally, the instantaneous helical vector) introduced in 
Eq. (2.401 allows one to vary the binary’s initial radial 


velocity, in addition to the control on the initial tangen¬ 
tial velocity provided by the orbital frequency D. If we do 
not make use of this freedom, i.e., if we set the radial ve¬ 
locity and the eccentricity parameter to zero, Vr = e = 0, 
we obtain the well known limit of standard quasicircular 
initial data. However, when we evolve such data, the sep¬ 
aration between the stars oscillates. From these oscilla¬ 
tions (or their imprint on the gravitational wave signal) 
we deduce actual measured eccentricities of e ~ 10“^. 
For the vast majority of astrophysical scenarios, this re¬ 
maining eccentricity is purely artificial, because GWs ef¬ 
ficiently circularise the orbit during the inspiral, leading 
to almost vanishing eccentricities during the last minutes 
before merger. In particular, the remaining eccentric¬ 
ity at merger e of the six observed double neutron star 
systems that will coalesce within a Hubble time will be 
e < 10-5 HI and thus several orders of magnitude smaller 
than the eccentricity obtained when evolving standard 
quasicircular initial data. 

The effect of this artificial eccentricity can be observed 
in various quantities, most notably the gravitational 
waveform. Therefore, we apply an iterative method to 
reduce the eccentricity, as outlined in m- This method 
is similar to the standard eccentricity reduction procedure 
for binary black holes HSHTHHZn], and the recent work 
on eccentricity reduction for binary neutron stars |29j . 


The basic idea is to find a measure for the eccentricity 
and determine which corrections have to be applied in 
order to remove the measured eccentricity from a Kep- 
lerian orbit. In this work, we use the proper distance d 
inside the hypersurface (measured along the coordinate 
line connecting the two local minima of the lapse, corre¬ 
sponding to the centers of the two stars) as well as the 
GW frequency UJ 22 to estimate the remaining eccentric¬ 
ity e. The coordinate distance dcoord can also be used to 
estimate e. However, dcoord depends much more strongly 
on the particular gauge choice than the proper distance 
does and thus gives reliable results only in certain cases 
(cf. |in])j which is why we choose not to use it. Using d for 
the two neutron stars, we track this quantity throughout 
the evolution and fit it to the model 

1 B 

d(t) = Agt — A.\t^ — — cos(cdft (j)) (4-9) 

2 (jjf 


Note that it is also possible to fit the time derivative d{t) 
and get rid of one fitting parameter. However, taking 
the derivative of d introduces noise, especially for lower 
eccentricities. Applying some low-pass filters can help to 
improve the results, but d is generally better suited for 
obtaining the eccentricity measure than d. The following 
results, however, apply to both methods. A comparison 
of the fitting model with the expected Keplerian orbits 
under the assumption that the eccentricity e is small for 
quasicircular orbits (so we can neglect higher orders of a 
series expansion), yields an eccentricity 


B 

doUJi 


■ 


(4.10) 


where c?o is the initial separation. 

The ellipticity part of the model ^ cos(a;ft-!-()))] leads 
to an initial radial velocity of B sin (p for the st ars, w hich 
means that the initial radial velocity in Eq. (2.40) has 
to be corrected by 


Svr = —B sin ( 


(4.11) 


in order to remove the eccentricity generated by the ra¬ 
dial velocity. Similarly, we can adjust the eccentricity 
parameter e (and thus the initial tangential velocity) in 
order to remove the residual eccentricity induced by ra¬ 
dial acceleration. The orbital angular frequency at apoas- 
tron Dapo is related to the D of the symmetry vector in 


Eq. (|2.40|) by Oapo « (1 — e)D. In Newtonian physics 
we have Dapo = (1 — e)GM/(ig. So if we change e 
by a small Se, this will lead to a change in Dapo by 
GM <■ 




-^6e = 


^6e. Here all contributions of 

e— 1 


O(e^) have been dropped, since we are only considering 
situations where is quite small (< IQ-"^). 

The radial acceleration in Eq. ( |4.9| ) is given by 
Bujf cos (p. Considering that a change of dUapo yields a 
change of the acceleration of do<5^apo! obtain a neces¬ 
sary correction of 


Bwfcosp 
02 n 

^apo^O 


(4.12) 
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where Hapo can be approximated by either ujf or Q,. We 


obtain the correction parameters (4.11), (4.12) from our 


fit to the data from the evolution, where we need at least 
one to two orbits for the fit to be accurate enough. Then 
we iterate this process until the eccentricity is sufficiently 
small, i.e., in most cases two or three iterations. 

Besides the already mentioned methods, there exist 
several other techniques to determine the eccentricity. To 
check the reliability of our results, we want to make use 
of the GWs to give an additional estimate of the eccen¬ 
tricity cgw- We follow the procedure given in [80], based 
on |5T] , and model the GW frequency motivated by post- 
Newtonian calculations as (cf. [55]) 


Wflt 




where ci, C 2 are determined by fitting and 


25M2 


T•o^ 


(4.13) 


(4.14) 


where tc and tq are again fitting parameters. We proceed 
by looking at the quantity 


eu,{t) = 


UJ22{t) - UJfjtjt) 


(4.15) 


This measure of eccentricity is time dependent and 
strongly oscillatory, but the global extremum 


CGw = niax I Ci 




(4.16) 


can be seen as a time-independent measure of eccentric¬ 
ity; see Fig. [^ As it is to be expected intuitively, this 
value can mostly be found in the beginning of the evolu¬ 
tion, just after the initial gauge noise has decayed. Note 
that the initial noise has to be cut off in order to obtain 
consistent data; this is also necessary for the estimate 
based on the proper distance. Gomparing the eccentric¬ 
ities obtained this way with the eccentricities computed 
from the proper distance, we observe an agreement within 
roughly five percent (comparing cgw.s = 8.4 x 10“'* to 
ed,3 = 8.7 X 10““* for the third eccentricity reduction iter¬ 
ation step of an SLy EOS run; cf. Table H- Considering 
the decrease in signal-to-(numerical)-noise of the proper 
distance oscillations at higher iterations, and thus the re¬ 
duced accuracy in the parameters one obtains from the 
fitting procedure, this agreement is satisfying and shows 
the consistency of the measures. In order to minimize 
the computational effort for finding eccentricity reduced 
initial data, it makes sense to evaluate to find the next 
iteration step’s parameters, since in this case we have to 
evolve for fewer time steps. Using the gravitational wave 
signal necessitates longer evolutions to ensure that the 
wave has reached the extraction radius. 

Table |TT| shows the numerical values for the eccentricity 
reduction iteration for two runs with different setups, viz., 
two equal-mass binaries with different total masses, one 
with the SLy EOS and the other with the r2 EOS. Addi¬ 
tionally, we show the proper distance d for the latter setup 
in Fig. We used SGRID with ua = ns = 2Q, = 8, 

ncart = 22 points. The evolution was done in BAM with 



FIG. 8: The eccentricity as a function of retarded time u 
computed from the residual of the GW frequency for SLy 
initial data with an initial [(2, 2) mode] GW frequency of 
Muj 22 ~ 0.0365. We shown the quasicircular data (red) and 
the third iteration of the eccentricity reduction (blue). The 
global extremum, marked by the points, gives eccentricities 
of eow.o = 0.0156 and ecw.s = 0.00084 for these interations, 
which is a factor of 20 improvement. The horizontal dashed 
lines mark the eccentricities Cd calculated from the proper 
distance. 


TABLE II: Iteration procedure for two binary setups starting 
at Muj 22 = 0.03 65. The first configuration is a r2 binary with 
individual masses of = 1.515 and a total ADM 

mass of Madm = 3.006. The second one is an equal mass SLy 
configuration with = 1.350 and Madm = 2.6782. 

The columns give the eccentricity e and radial velocity Vr 
input to the code [cf. Eq. ( |2.40[ )], along with id, the remaining 
eccentricity in the binary evolution measured using the proper 
distance d in the hypersurface. We also list the values of the 
binding energy Eb = Madm — M and the angular momentum 
Jadm, which we normalize by M and M^, respectively. 


EOS 

Iter 

e [10-3] 

Vr [10“3] 

id [10-3] 

Eb/M [10-3] 



0 

0 

0 

9.77 

-7.984 

1.0700 

r2 

1 

-6.8 

-1.63 

1.38 

-7.922 

1.0729 

2 

-5.7 

-1.14 

0.91 

-7.920 

1.0738 


3 

-6.3 

-1.16 

0.56 

-7.920 

1.0734 


0 

0 

0 

12.41 

-8.115 

1.0541 

SLy 

1 

-6.0 

-1.13 

7.80 

-8.103 

1.0580 

2 

-12.1 

-1.91 

3.97 

-8.088 

1.0615 


3 

-13.7 

-1.09 

0.87 

-8.085 

1.0625 


a constraint damping Z4c evolution scheme, as described 
in Sec. |V| We used a total of 7 refinement levels, where 
the two inner levels contained 96 points in each direction 
in each box with a grid spacing of 0.15 in the finest. The 
other outer boxes all contain 192 points in each direc¬ 
tion; the grid spacing is doubled on each successive level 
moving outwards. Finally, the outermost level is given by 
a cubed sphere (cf. jM] [55H55] ; this is also used for the 
evolutions in Sec. with 192 and 84 points in the radial 
and azimuthal directions, respectively. Over the course 
of three iteration steps, using the method outlined above 
to iteratively determine the correction parameters and 
Vr, we were able to decrease the eccentricity in the r2 case 
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FIG. 9: Comparison of the eccentricity of the simple poly¬ 
tropic r2 setup by looking at the proper distance as a function 
of time. The blue dotted line is from an evolution of the orig¬ 
inal initial data set, while the solid black line uses the third 
iteration step of the eccentricity reduction procedure. 


from ed.o = 9.8 x 10“^ to ed, 3 . = 5.6 x 10“^, which can 
be seen in a significant improvement of the oscillations 
of the proper distance (Fig. |^. Similar improvements in 
the gravitational waves can be seen in Fig. for the SLy 
case. 

Furthermore, we can try to give some PN estimates of 
the improvement of the data by looking again at the PN 
expressions from Mora and Will |76] . which we already 
utilized in Sec. |IVB| We can even expand our compar¬ 
isons to fourth post-Newtonian order (as summarized in 
[29]) without any essential changes in the results, since 
the 4PN contribution is small. 

Given the fairly large initial separation of the bina¬ 
ries, we expect the fourth post-Newtonian order results 
for the ADM energy and angular momentum in terms 
of the binary’s angular velocity to be quite accurate and 
we indeed find that the eccentricity reduced data gives 
a better match to the values for these quantities for a 
circular orbit than do the original data: Table |n| gives 
the binding energy and the angular momentum Jadm 
and if we compare these to PN for the SLy run (with 
Jadm, pn = 1.0637Af^ and E^^ pn = —0.008081Af), we 
find a relative error of ^ 1% for original data and an error 
of ~ 0.2% for the data at the third iteration. Similarly, 
for the simple polytropic r2 run, we compare our results 
to Tadm, pn — 1.0731iW^ and E^^ pn = —0.007918117 and 
find an improved agreement of the eccentricity reduced 
data with a deviation of ^ 0.04%, compared with a devia¬ 
tion of ~ 0.4% for the original data. We thus see that the 
eccentricity reduced data are closer to quasi-equilibrium 
than the standard initial data, as was noted by |29j . 


D. Unequal masses 

As outlined in Sec. and explained in more detail in 
Appendix |A 1[ some population synthesis models pre¬ 
dict relatively high mass-ratio BNS systems. In this 
section we want to illustrate SGRID’s ability to gener¬ 
ate quasiequilibrium configurations for such systems. We 
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FIG. 10: q — 2.06 sequence. Top panel: ADM mass compared 
with 4PN tidal results (dot dashed lines). Bottom panel: 
ADM angular momentum compared with 4PN tidal results 
(dot dashed lines). 


model the systems by the stiff EOS MSlb. This EOS 
allows rather large baryonic masses up to Mf, = 3.35 and 
gravitational masses up to M = 2.76 for single neutron 
stars in isolation. As is seen in the next subsection, it is 
one of the EOSs for which we achieved high neutron star 
masses in equal-mass binaries with SGRID. 


The particular configuration constructed below con¬ 
sists of = 1.00 (M^ = 0.94) and = 2.20 (M^ = 
1.94) neutron stars resulting in a mass ratio oi q = 2.06. 
We want to highlight that this is (i) a slightly larger mass 
ratio than the largest predicted b y th e population synthe¬ 
sis models discussed in Appendix ] A 1[ (ii) the largest mass 
ratio computed for a realistic (irrotational) binary config¬ 
uration in quasi-equilibrium,® and (iii) the lar gest m ass 
ratio evolved in full general relativity;® see Sec. V A We 
employ a grid with ua = tib = 28, n^p = 8, ncart = 24 
and construct a sequence varying the distance param¬ 
eter b G [16,30]. This results in GW frequencies of 
Moj22 G [0.019,0.041]. 


The ADM mass and angular momentum for this se¬ 
quence are shown in Fig. jlOj As a comparison, we show 
the 4PN results including tidal components given in Ap¬ 
pendix A of [H] (obtained from |H5] and ^Tj). However, 
the influence of the tidal contributions and higher-PN 
terms (above 2PN) is negligible at the scale shown here. 
Additionally, we emphasize that a direct comparison of 
these results is not really warranted, since SGRID uses 


® In |22| (see that paper’s Table 1), one of us computed corotating 
binary neutron star initial data with baryonic mass ratios up to 3 
for large separations, using a polytropic EOS and small baryonic 
masses. 

® Up to now the highest mass ratio evolved in full general relativity 
was q = 1.5 |35II36| . 
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the conformal flatness approximation, which is known to 
be violated at 2PN (see, e.g., the discussion in Sec. Ill A 
of [17]). 

Due to the newly implemented mechanism in SGRID 
which adjusts the center of mass so that T^dm i® 
small (ideally vanishing), the linear momentum of the 
configuration stays within ~ 10“^ to ^ 10“^. In par¬ 
ticular, for the configuration with Mlo22 — 0.03 59 that 
we evolve in Sec |V A| the linear momentum is Padm = 
( 4.23 X 10 "'^, - 1.10 X 10 - 5 ,- 1.96 x 10 "®). 


E. High compactness 


While population synthesis models predict a num¬ 
ber of binary neutron star systems containing high- 
mass ( and t hus high-compactness) neutron stars (see Ap¬ 
pendix |X^, computing initial data for binary neutron 
stars with high compactnesses is a challenging task for 
most codes. The highest compactness achieved for a neu¬ 
tron star in a relativistic binary is: C ~ 0.26 for a neutron 
star-black hole system jSS], C ~ 0.25 for corotating binary 
neutron stars [89] . C ~ 0.26 for irrotational binary neu¬ 
tron stars m, and C ~ 0.22 for binary neutron stars with 
relatively low spins [18] . 

In the following we present results for most EOSs listed 
in Table H] (except for the r2 EOS, which does not allow 
high compactness neutron stars) and provide an estimate 
of the highest compactness easily reachable with our nu¬ 
merical method. In physical terms, the maximum feasible 
compactness depends mostly on two different properties: 
(i) the chosen EOS and (ii) the binary separation. How¬ 
ever, in SGRID the resolution and iteration procedure 
(e.g., the softening parameter and mass increment) also 
play important roles. 

As a starting point, we choose a simple procedure to 
estimate the maximum compactness C„iax reachable with 
SGRID. We consider equal-mass nonspinning binaries, 
fix the separation parameter 6, and increase the bary- 
onic masses of the two stars in steps of AMf, = 0.005, 
starting from Mj^ = = 1.2. We stop increasing the 

mass when the elliptic solve of Eqs. (2.18), in particular 


the Hamiltonian constraint equation, does not reduce the 
residual and no solution can be found. For our test we 
use a resolution of ua = ns = 24, = 8, ncart = 20. 

Table jlllj summarizes the maximum compactness 
achieved for different EOSs and separations. In all cases 
we were able to achieve higher compactnesses at larger 
separations. The maximum compactness for the EOSs 
and separations considered lies within (0.197,0.212) for 
this simple iteration procedure. We obtain the highest 
compactness for this test (with a constant number of 
grid points) with the simple polytropic EOS r2.72. This 
is to be expected, since runs with piecewise polytropic 
EOSs require higher resolutions to obtain the same accu¬ 
racy as runs with simple polytropic EOSs, as discussed 
in Sec. IIVFI 

Even higher compactnesses can be achieved by reduc¬ 
ing AMb, decreasing the softening parameter to 0.2, and 
using higher resolutions. We now give an explicit example 
of this procedure. 


In order to achieve high compactness with the SLy 
EOS, we have considered a non-spinning equal mass sys¬ 
tem with a separation parameter of 6 = 31. We started 
with baryonic masses = 1.5 and increased 

the mass in each iteration by a factor of 1.07 until we 
reached = 1.82. The highest resolution used 

for this configuration was ua = ns = 28, Utp = 8, 
^^Cart = 24. This results in a binary with ADM mass 
Madm = 3.193, angular momentum Jadm = 12.759, or¬ 
bital angular velocity 11 = 0.00336, dimensionless GW- 
frequency Muj 22 = 0.02 1 6, and a coordinate separa¬ 
tion of 63.270 between the star centers. A single TOV 
star with the same baryonic mass would have an ADM 
mass of 1.60564 and a radius of 7.69608 (in standard 
Schwarzschild areal radius coordinates), which implies a 
compactness of 0.2086. 

Increasing the baryonic masses beyond = 

1.82 has proven very difficult. For masses above this value 
we observe that the elliptic solver cannot always solve 
Eq. (2.18a) for ip within the main iteration, because the 
values from the previous iteration are not good enough 
as an initial guess for the Newton-Raphson scheme we 
use. We were able to address this problem by making 
two changes to our Newton-Raphson scheme. First, we 
take smaller Newton steps if this results in a smaller 
residual error than a full step (this procedure is known 
as backtracking). And second, if backtracking also fails 
we simply skip the elliptic solve for pi. The overall it¬ 
eration still succeeds if we do not skip the elliptic solve 
for p too often. Using this trick we started from the 
initial data for = 1.82 and increased the 

masses by a factor of 1.02 in the main iteration un¬ 
til we reached = 2.0. The resulting ini¬ 

tial data were again computed with ua = tib = 28, 
= 8, ncart = 24 and result in a binary with ADM mass 
ATadm = 3.459, angular momentum Jadm = 14.444, or¬ 
bital angular velocity D = 0.00347, dimensionless GW- 
frequency Muj 22 = 0.0242, and a coordinate separation 
of 63.262 between the star centers. The Hamiltonian con¬ 
straint violation for this data set is about twice as big as 
for the lower mass data, because we do not always solve 
Eq. (2.18aI, which is the Hamiltonian constraint. A sin¬ 
gle TOV star with the same baryonic mass would have an 
ADM mass of 1.74067 and a radius of 7.6104 (in standard 
Schwarzschild area radial coordinates), which implies a 
compactness of 0.2287. 

If we try to increase the baryonic masses beyond = 
= 2.0 we can never solve Eq. ( [2.18a] ) for p, so we 
obtain constraint violating initial data. 


F. Convergence of SGRID 

The convergence of SGRID was already presented in 
Fig. 3 of [22] for a corotating equal-mass quasicircu¬ 
lar binary neutron star system with a simple polytropic 
EOS. However, we want to show here how well the 
code converges in more complicated situations. To il¬ 
lustrate this, we consider the configurations presented 
in Tab. m In particular, we choose configurations with 
(7 = 1.0 where each star has a gravitational mass of 
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TABLE III: Overview of high compactness conhgurations without fine tuning the iteration procedure. The columns refer to: 
the SGRID distance parameter b, the dimensionless gravitational wave frequency Muj 22 , the maximum baryonic mass Mb, max, 
the corresponding gravitational mass Mmax, and compactness Cmax- 


EOS 

b 

MuJ22 

Mb, max 

-^^max 

^max 

b 

MlJ22 

Mb, max 

Mmax 

^^max 

b 

MuJ22 

max 

Mmax 

Cmax 

SLy 

22 

0.032 

1.727 

1.534 

0.199 

26 

0.026 

1.744 

1.547 

0.200 

32 

0.020 

1.780 

1.575 

0.204 

ENG 

22 

0.033 

1.788 

1.583 

0.199 

26 

0.024 

1.824 

1.611 

0.203 

32 

0.021 

1.852 

1.632 

0.206 

MPAl 

22 

0.035 

1.871 

1.661 

0.200 

26 

0.029 

1.899 

1.676 

0.202 

32 

0.022 

1.927 

1.683 

0.203 

ALF2 

22 

0.036 

1.899 

1.678 

0.200 

26 

0.029 

1.937 

1.708 

0.204 

32 

0.023 

1.976 

1.738 

0.208 

H4 

22 

0.038 

1.976 

1.762 

0.197 

26 

0.031 

1.996 

1.777 

0.199 

32 

0.024 

2.026 

1.800 

0.203 

MSlb 

22 

0.043 

2.183 

1.931 

0.198 

26 

0.036 

2.216 

1.963 

0.201 

32 

0.027 

2.250 

1.982 

0.203 

MSI 

22 

0.043 

2.205 

1.955 

0.197 

26 

0.036 

2.238 

1.981 

0.200 

32 

0.028 

2.273 

2.007 

0.203 

r2.72 

22 

0.042 

2.162 

1.908 

0.205 

26 

0.034 

2.183 

1.924 

0.207 

32 

0.027 

2.238 

1.966 

0.212 


TABLE IV: Convergence test setups. Columns refer to: 
Name of the configuration, EOS, distance parameter b, ec¬ 


centricity parameter e [from (2.38l], the baryonic masses 


and Mb , and the angular velocity vector a;* as specified in 


Eq. (2.221. 

Name 

EOS 

6 

e 

Mb^ 

Mf 

<^A,B 

r2ql00w000e00 

r2 

20 

0.0 

1.4336 

1.4336 

(0,0,0) 

r2ql00w005e00 

r2 

20 

0.0 

1.4336 

1.4336 

0.005- (1,1,1) 

r2ql00w000e03 

r2 

20 

0.3 

1.4336 

1.4336 

(0,0,0) 

r2ql00w005e03 

r2 

20 

0.3 

1.4336 

1.4336 

0.005- (1,1,1) 

r2qll6w000e00 

r2 

20 

0.0 

1.5491 

1.3199 

(0,0,0) 

H4ql00w000e00 

H4 

20 

0.0 

1.4687 

1.4687 

(0,0,0) 

H4ql00w005e00 

H4 

20 

0.0 

1.4687 

1.4687 

0.005- (1,1,1) 

H4ql00w000e03 

H4 

20 

0.3 

1.4687 

1.4687 

(0,0,0) 

H4ql00w005e03 

H4 

20 

0.3 

1.4687 

1.4687 

0.005- (1,1,1) 

H4qll6w000e00 

H4 

20 

0. 

1.5887 

1.3506 

(0,0,0) 


1.35. We either investigate nonspinning configurations 
[i.e., uj\ = ujg = (0,0,0)] or choose an angular velocity 
of uj\= ojg = (0.005,0.005,0.005). We also use different 
eccentricities of e = 0.0 and e = 0.3. Finally, we consider 
an unequal mass configuration with gravitational masses 
1.45 and 1.25 {q = 1.16). We compute initial data for 
these systems for both a simple polytropic EOS as well 
as a piecewise polytropic EOS (the fit to the realistic H4 
EOS). 

Figure [TT] summarizes our findings. In the simple poly¬ 
tropic runs (upper panel), we encounter a higher Hamilto¬ 
nian constraint for r2ql00w005e00 and r2ql00w005e03 
than for the other configurations. This clearly shows that 
the addition of spin increases the complexity of the sys¬ 
tem, while adding eccentricity has no noteworthy effect. 
In the cases where no spin is present, the Hamiltonian 
constraint is one to two orders of magnitude smaller. For 
unequal masses we see a slightly larger constraint vio¬ 
lation. The momentum constraints converge in a similar 
manner, but their magnitude is roughly one order of mag¬ 
nitude smaller. 

When piecewise polytropes are employed, the con¬ 
straint violations increase substantially. This is to be 
expected, since the solution is only at the interfaces 
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FIG. 11: Convergence analysis for representative configura¬ 
tions. Simple polytropic r2 configurations (top panel) and 
piecewise polytropic (H4) configurations (bottom panel). The 
plot shows the average Hamiltonian constraint inside the two 
stars in the regions A € [0, Amax] and thus includes the stars’ 
surfaces, which are the most problematic regions. We have 
fixed = 8 and use ua = ub, ncart = ua — 4. 


between the different pieces, so one will no longer obtain 
exponential convergence from the spectral method. In 
this case, however, the addition of spin and eccentricity 
or the consideration of unequal masses has no noteworthy 
additional effect. Since spin most likely affects the con¬ 
vergence primarily through the deformation of the star’s 
surface, it is not surprising that it does not affect the con¬ 
vergence so much for piecewise polytropes at these resolu¬ 
tions, where the portion near the surface is already quite 
troublesome, due to the necessity of resolving the crust. 
We expect that one would see a difference due to the 
addition of spin at much higher resolutions or for higher 
spin magnitudes. This reduction in convergence when us¬ 
ing piecewise polytropes means that in most cases we use 
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higher resolutions when computing piecewise polytropic 
setups, compared to previous work with just simple poly¬ 
tropes. 


V. DYNAMICAL EVOLUTIONS 

In this section we present evolutions of our new 
initial data for two new configurations: (i) a high-mass 
ratio simulation with q = 2.06, which is the highest 
mass ratio binary neutron star ever evolved in full 
numerical relativity; (ii) an unequal-mass configuration 
where the spins are not aligned with the orbital angular 
momentum, which is the first processing binary neutron 
star simulation performed so far. The new dynamical 
evolutions presented here show, as a proof of principle, 
that the BNS phenomenology can be very rich in these 
newly accessible regions of parameter space, where we 
are now able to study mass exchange during the inspiral 
or the precession and nutation of the orbital plane 
during a binary neutron star simulation. Additionally, 
we discuss (iii) the effect of eccentricity reduction on the 
waveform phasing. 

We perform the evolutions with the newest version of 
the BAM code [36l EB ES] including the recently im¬ 
plemented conservative mesh refinement algorithm em¬ 
ployed in EH ED]. We also use the Z4c scheme El ED 
with constraint preserving boundary conditions El E2] ■ 
The BAM grid setup consists of 7 refinement levels. The 
outermost level {I = 0) uses a multipatch “cubed-sphere” 
grid EIPPP). 

The setups of Sec. lYl do not employ any symmetry 
condition. For the other simulations we employ reflection 
symmetry across the orbital plane, letting us use half the 
number of grid points in the z-direction. Further infor¬ 
mation about the particular grid and initial conditions is 
given in Tab. jVj 


A. A q — 2.06 binary 



-2.500 

-3.379 

-4.259 

-5.138 

-6.017 


-8.655 

-9.534 
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FIG. 12: Plot of the matter density po and the velocity a* 
in the orbital plane for the q = 2.06 simulation at times t = 
1726M (upper panel), two revolutions later at t = 2227M 
(middle panel), and yet another two revolutions later at t = 
2644M (lower panel). Over the 4 revolutions shown, one can 
see a clear mass transfer between the two stars. Note that 
each plot has a somewhat different scale. 


In this section we evolve one of the models computed 
in Sec. |IVD| We have set the rest mass of the primary 
star to Mff = 2.200, corresponding to a gravitational 
mass of = 1.944 and compactness = 0.199 
in isolation, while the companion is characterized by 
= 1.000, = 0.944, C® = 0.103. We have not 

tried to add spin/eccentricity or performed eccentricity 
reduction since we want to focus solely on the high mass 
ratio effects. Nevertheless, the setup shows a rather small 
eccentricity of = 4.2 x 10“^, cgw = 5.1 x 10“^ and 
despite the high mass ratio a small linear momentum of 
Padm = (4.23 X 10-^ -1.10 X -1.96 x lO"®). The 
initial gravitational wave frequency is MUJ 22 = 0.0359. 
Star B is already deformed at this separation indicated 
by a mass shedding parameter of y ~ 0.89; cf. Eq. (C4). 

This configuration, with a gravitational mass ratio of 
q = 2.06, is the highest mass ratio binary neutron star 
evolved in full general relativity, and well above the 
mass ratio of q = 1.5 considered before EB ESj. Note 


that E3] reports results of Newtonian smooth particle 
hydrodynamics (SPH) evolutions (including radiation re¬ 
action) for a similar q = 2 setup [gravitational masses 
of (1.0 -I- 2.0)], but does not mention any mass transfer, 
which we find in our simulation. However, mass transfer 
is observed in white dwarf binary simulations with q = 2 

m- 


1 . Mass transfer 

In Fig. we show snapshots of the density in the 
orbital plane during the inspiral. The upper panel of 
Fig. shows the binary at t = 1726M. Although the 
stars are clearly separated, mass transfer from the com¬ 
panion (M^) to the primary star (M^) has already set 
in. Two revolutions later [t = 2227M (middle panel)] 
and four revolutions later [t = 2644M (lower panel)] the 
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TABLE V: Initial data and grid details for the dynamical evolutions. The columns refer to: the simulation name, the 
gravitational and baryonic masses of star A and B, the stars’ dimensionless angular momenta , the number of grid points 
employed in each direction in the fixed and moving boxes, the number of radial and azimuthal points, the Hnest grid spacing, 
and the outer boundary location. 


name 






3^ 

n 

nrav 

Ur 

ne 

/ig 

Tb 

MSlb-q206 

1.944 

2.200 

0.944 

1.000 

(0,0,0) 

(0,0,0) 

128 

72 

144 

63 

0.250 

1692 

SLy(^^) 

1.3553 

1.500 1.1072 1.200 (0.13/y3)(l,l,l) (0.16/y3)(l, 1,1) 

128 

64 

128 

56 

0.245 

1532 

SLy(^^) 

1.3547 

1.500 

1.1067 

1.200 

0.077(0,0,1) 

0.089(0,0,1) 

128 

64 

128 

56 

0.245 

1532 

SLy(°°) 

1.3544 

1.500 

1.1065 

1.200 

(0,0,0) 

(0,0,0) 

128 

64 

128 

56 

0.245 

1532 

SLy-eccred 

1.350 

1.495 

1.350 

1.495 

(0,0,0) 

(0,0,0) 

192 

96 

192 

84 

0.15 

1384 


mass transfer becomes more dramatic. During this pe¬ 
riod, ~ (2-3) X 1O“^M0 of material was transferred be¬ 
tween the two stars, i.e., ^ (2-3)% of the rest mass of 
the less massive star. The estimate is based on the rest 
mass leaving the finest refinement level around star B and 
entering the refinement level of star A. The uncertainty 
is mainly related to mass loss due to the artificial atmo¬ 
sphere treatment [35] ■ (Note that the overall rest mass 
is conserved to better than 0.12% in this simulation until 
very late times, post-merger, when matter starts leaving 
the grid.) We observe that for this system the mass trans¬ 
fer happens continuously until the companion is tidally 
disrupted. 

The average rate of mass transfer is Mab 10“® ~ 
IMq s“^, taking place for ^ 10“^ s, from which one can 
estimate the accretion power. Here we just make a sim¬ 
ple, order-of-magnitude estimate, since this is likely all 
that is warranted by the accuracy of our estimate of the 
mass transfer. If one just considers the change in en¬ 
ergy of this matter going down the more massive star’s 
Newtonian potential well, we have an average accretion 
power of ^ MabC^ ~ ~ 10®^ erg s“^, compara¬ 

ble to the neutrino luminosities found in simulations of 
BNS mergers [HI |49l |l0[. (Recall that is the com¬ 
pactness of in isolation.) This gives a total accretion 
energy of ~ 10®^ ergs, which is comparable to the energy 
released in a supernova. But of course, one cannot say 
anything definite about the amount of this energy that 
would be released in photons or neutrinos, since these are 
not present in our simulation. While one might expect to 
be able to see the effects of this amount of mass transfer 
on the phase of the gravitational waves, the present sim¬ 
ulation does not appear to be accurate enough to make 
such a comparison. 

The merger happens at f = 2692M. We classify the 
merger remnant as a supramassive neutron star (SMNS), 
since it is below the maximum supported gravitational 
mass of a rigidly rotating star for the MSlb EOS, which is 
^ 3.2-3.3M0, i.e., roughly ^ 15-20% larger than the cor¬ 
responding maximum gravitational mass of a nonrotating 
neutron star (see, e.g.. Sec. 2.9.1 in [95]). We thus do not 
expect the SMNS to collapse on dynamical timescales. 
This is supported by our simulation, where no indication 
for a collapse is present through the end of the simula¬ 
tion at < = 6500M ~ 0.09 s. The central density reaches 
a constant value — 9.8 x 10“^ ~ 6 x 10^^ g cm“^ and 
the final merger product settles into a stable configura- 
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FIG. 13: Plot of the matter density po, the density of unbound 
matter pou, and the velocity field u® in the orbital plane for 
the q = 2.06 simulation after merger. In the upper panel 
(t = 2798M) a clear spiral like pattern in the ejecta is visible, 
where material is expelled due to torque in the tidal tail of 
the companion star. The lower panel (342IM) shows that 
material is ejected over the entire grid anisotropically. Note 
that the two panels have different scales. 


tion. 


2. Ejecta and kick 

In this simulation we observe a significant mass ejec¬ 
tion, Mejecta — 7.6 X Mq, which is among the largest 
found for full general relativistic simulations of binary 
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neutron star mergers, including the case of eccentric bi¬ 
naries |96] . and much larger than any of the ejecta masses 
found for the quasicircular case in the studies in [36l EZ] ■ 

Figure[^ visualizes the ejected material, distinguishing 
it from the bound material by using a different color bar. 
Here we compute the unbound material using the method 
given in [36] . In our simulation most of the material is 
ejected into the orbital plane by torque on the tidal tail of 
the companion star; cf. the discussion in [SiET]. We can 
also see this in the spiral-like pattern in the upper panel of 
Fig. at a time of t = 2798M. In the lower pattern we 
see also clearly that the ejection happens anisotropically, 
where the density inside of the ejected material at a given 
radial distance from the SMNS differs by several orders 
of magnitude in different directions. The kinetic energy 
of the ejecta is ^ 2 x 10“^ ~ 4 x 10®° erg. Note that 
the ejecta mass we find is a factor of ^ 2 larger than 
that found in the Newtonian calculation of [93] (Table 1 
there). This is significant, since, as discussed in [36] . 
the uncertainty on the ejecta mass is about < 20%, and 
mainly due to the resolution. 

The anisotropic mass ejection causes the merger rem¬ 
nant to recoil. We approximate the ejecta’s linear mo¬ 
mentum by computing 

Pei = Mej (-Uplane) = Mej J B^dy ’ 

where the integrals here are restricted to the orbital 
plane (while Mgj is computed from the integral of the 
rest-mass density D of the unbound matter over all three 
dimensions) and Upiane denotes the ejecta velocity in the 
orbital plane (i.e., just the x and y components). We find 
that = ||Pej||/M ~ 100—1000 km s“^. This number 
is also consistent with the one obtained from the coordi¬ 
nate position of the SMNS. However, the value should 
just be seen as an order of magnitude estimate, where 
the main difficulties here are the low resolution, the long 
simulation time, and the gauge dependence of the mea¬ 
surement. Since this is an unequal-mass system, we also 
expect some contribution to the kick from anisotropic 
gravitational wave emission. We compute this from the 
GW linear momentum flux, 

=-^ J Powdt (5.2) 

where the linear momentum flux Pqw is computed as in 
|KT] . We find ~ 100 km s“^ at merger. The GW 
kick is smaller than the ejecta kick, as is found in black 
hole-neutron star simulations [551 EH] ■ 



FIG. 14: The four dominant multipoles of the GWs from the 
q = 2.06 simulation, viz., (l,m) — (2,2), (2,1), (3,3), and 
(4, 4). We plot the real part of the modes in black and show 
the dimensionless GW frequency for all modes in red. In the 
top plot, the solid vertical line marks the merger. The fre¬ 
quency oscillations are in part unphysical and resolution de¬ 
pendent. The large frequency spikes in the post-merger phase 
are caused by zeros of the amplitude. 



/[kHz] 

FIG. 15: Power spectrum density for the dominant postmerger 
modes of the q = 2.06 simulation. We used a time interval u £ 
[3800M, 5700M] for the computation, and therefore focus only 
on the SMNS spectrum. The harmonicity of the frequencies 
fi, f 2 , fa, fi is clearly visible. 


3. Gravitational waves 


Let us discuss the gravitational waveform. The four 
dominant modes of the gravitational waveform are pre¬ 


sented in Fig. 14 The inspiral-merger signal ends at 
Mmrg = 2692M, and we find no evidence for an obvious 
GW signature of the mass transfer described above. We 
study the effect of the mass-ratio on the GW multipolar 


structure by computing the relative GW energy contribu¬ 
tion of each dominant mode over the total energy released 
up to merger {Eim/E) and comparing with the q = 1 case 
with the same EOS [l(K)j (though with a slightly different 
total mass: 2.70, as opposed to 2.89 for the q = 2.06 sim¬ 
ulation; the q = 1 system also radiates almost 40% more 
energy per total mass than the q = 2.06 system). We find 


A/W; 
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that for the q = 1 simulation ^ 99.6% of the energy is 
released in the (2, 2)-mode, while ^ 98% is released for 
MSlb-q206. For q = 1 the (2,1) and the (3, 3) modes are 
zero by symmetry; for q = 2.06 the ( 2 , 1 ) mode releases 
about ^ 0.08% of the energy and the (3,3) mode ^ 1.4%. 
The next dominant mode, (4,4), contributes to ~ 0.19% 
of the total energy for q = 1 and 0 . 22 % for q = 2.06. 

Finally, we study the postmerger GW frequencies for 
the dominant modes. We perform a Fourier analysis on 
the interval u S [3800M, 5700M] in Fig. 15 The /2 fre¬ 
quency [the dominant frequency of the ( 2 , 2 ) mode] is 
clearly visible in all modes (though it may not be dom¬ 
inant) and is close to /2 = 2.09kHz (M /2 = 0.0298). /2 
agrees within 2.5% with the relation presented in [lOlj . 
which indicates the robustness of the fit given there even 
for high mass ratios. Comparing with the q = 1 case, 
we do not find a significant difference in the /2 value. 
Inspection of the peaks of the other modes, reveal that 
the fk frequencies are harmonic to a very high accu¬ 
racy, i.e., fi = ( 2/2 = fs/S = (4/4:. In particular we 
have /i = 1.05 kHz (M/i = 0.0150), h = 3.16 kHz 
(M /3 = 0.0449), /4 = 4.20 kHz (M /4 = 0.0598). The 
agreement is better than 1 %, though the uncertainties 
in the frequencies are larger than this, about 0.15 kHz 
and are mainly caused by a shift of the frequencies over 
time due to the changing compactness of the neutron 
star. This harmonicity was also found in |102) . though 
there they obtained the mode frequencies from a Fourier 
transform of the pressure and used the spatial conformal 
flatness approximation in their simulation. Here we have 
verified that the same harmonicity is present in an observ¬ 
able quantity (the modes of the gravitational radiation) 
in a fully relativistic simulation. 


B. A precessing unequal mass binary 

In this section we investigate the evolution of a BNS 
system with a rather generic spin configuration and q = 
1.22. The CRV angular velocity w® for star A and star B 
is set as uj = 0.005(1,1,1), i.e., their spins both point at 
an angle of 45° to the orbital angular momentum. This 
evolution is the first precessing BNS merger simulated in 
numerical relativity. We indicate this precessing spin con¬ 
figuration with SLy^'^^) as we use the SLy EOS (Tab.|^. 
The dimensionless spin magnitudes are ~ 0.13 and 
~ 0.16. Although such high spins are not observed in 
double neutron star systems so far, there is no physical 
reason to exclude such a scenario. In particular, binary 
neutron star systems with spins misaligned from the or¬ 
bital angular momentum, as here, are more likely to form 
dynamically in dense stellar regions, such as globular clus¬ 
ters, where there are many rapidly spinning neutron stars, 
as is discussed in Appendix |A 2| The stars have baryonic 
masses of = 1.5 and M/* = 1.2; the gravitational 
masses are M^ = 1.3553 and = 1.1072. 

Together with SLy^^'^^ we evolve for comparison two 
other configurations with the same baryonic masses (and 
q) but with the CRV angular velocity set to a; = 
0.005(0,0,1) (SLy%f)), and w = 0.0 (SLy(oo)). Thus 
SLy*^f^^ has no precession, and SLy*^°*^) no spin interac¬ 


tions at all. 

All the SGRID data are computed with ua = ns = 28, 
= 8 , ncart = 24. We have not tried to reduce the ec¬ 
centricity for this particular setup to save computational 
resources. The residual eccentricity of the precessing sys¬ 
tem is Cd ~ 4 X 10“^, CGW — 5 X 10“^ for the two eccen¬ 


tricity measures considered in Sec. IV C The resolution 
in the Hnest level covering the neutron star is hg = 0.245, 
similar to the low resolution setups of |100Lll03[ll04j . but 
no symmetries are applied to the grid (full 3D). Thus, al¬ 
though the principal dynamics are properly modeled at 
this resolution, quantitative statements come with quite 
large uncertainties. 


1. Dynamics 


Figure 16 (upper and lower left) shows the coordinate 
tracks (positions of the local minimum of the lapse) of the 
neutron stars in the SLy^^^^ simulation. The change of 
the orbital plane due to the misaligned spin of the binary 
neutron stars is clearly visible. This effect can also be 
seen in the change of the z-coordinate over time (upper 
right panel). During the inspiral approximately one pre¬ 
cession cycle is finished, i.e., the orbital plane again coin¬ 
cides (approximately) with the xy-plane during merger. 

We present the evolution of the norm of the Hamil¬ 
tonian constraint in the upper panel of Fig. [iT] and of 
the rest mass conservation in the lower panel. Because 
of the constraint propagation and damping properties of 
the Z4c scheme, the constraint violations stay at or below 
the value of the initial data. The rest mass conservation 
over the entire simulation is up to ~ 0.3%. Overall, these 
diagnostics indicate the errors in the simulation are under 
control. The small violation of rest-mass conservation is 
related to the artificial atmosphere and the relatively low 
resolution we employed; see the discussion in |36j . 

The initial angular momentum of the system is 
JADuit = 0) = (0.251,0.239,6.951) or normalized J = 
(0.0361,0.0343,0.9989). We find that J changes slightly 
over the entire simulation, which is to be expected, as 
the total angular momentum precesses slightly in post- 
Newtonian calculations, as discussed in, e.g., |1()5] . We 
can estimate the opening angle of the precession cone (ne¬ 
glecting spin-spin effects) using Eq. (51) from [105] and 
evaluating everything at < = 0 , using the initial orbital 
frequency to evaluate M/r to IPN order. We calculate 
the orbital angular momentum by subtracting the spin 
angular momenta of the two stars in isolation from the 
ADM angular momentum. Using this method, we obtain 
(as a hrst approximation) an opening angle of the pre¬ 
cession cone of Aj ~ 1.3 x 10“^. Due to the decreasing 
orbital separation over the evolution, we expect that the 
opening angle should be slightly larger in our full GR 
simulation. In fact, this can be observed and we find an 
angle of ^ 1.5 x 10“^ ~ 0.086°, where the numerical un¬ 
certainty is < 10 “^ based on comparison with simulations 
without precession. The opening angle of the precession 
cone for the orbital angular momentum is Al ^ 0.05, 
from a simple calculation of the initial angle between the 
orbital and total angular momenta, which agrees with the 
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FIG. 16: Orbital dynamics of the precessing binary neutron star inspiral. Upper left: neutron star tracks in the a;i/-plane. Upper 
right: ^-coordinates of the stars’ centers as a function of the coordinate time t. Lower left: 3D neutron star tracks visualizing 
the precession of the orbital plane, where the spheres denote the original positions of the stars. Lower right: Precession and 
nutation of the orbital angular momentum L (green lines). The coordinate system is rotated such that J(t = 0) lies along the 
z-axis. The orbital angular momentum performs slightly more than one precession cycle with a period of Tprecess — 4720M. The 
opening angle Al decreases from ~ 0.05 ~ 3° to ~ 0.035 ~ 2° over the inspiral. The direction of the total angular momentum 
also processes at the same period with a considerably smaller opening angle of Aj ~ 1.5 x 10“® ~ 0.086° (red lines). 


initial opening angle found in our simulation. The open¬ 
ing angle Xl decreases from ^ 0.05 ~ 3° to ^ 0.035 ~ 2° 
over the inspiral. 

The lower right panel of Fig. |16| presents the precession 
of the orbital plane, as given by the direction of the or¬ 
bital angular momentum. The effects of precession and 
nutation are clearly visible. The normalized orbital an¬ 
gular momentum L we plot in this figure is constructed 
as the vector orthogonal to the orbital plane, which we 
estimate using the coordinate line between the two star 
centers at two adjacent timesteps. To minimize high fre¬ 
quency noise we apply a low-pass filter. The precession 
period is Tprecess — 4720M, which agrees within < 10% 
with the PN estimates, based on a PN evolution similar 
to the discussion in IZO]." 


Binding energy vs. reduced orbital angular momentum 
curves provide a gauge-invariant way of characterizing 
the binaries’ dynamics [551 [Ml ITUni. In order to compute 


such curves, we modify Eqs. (4.1) and (4.2) to take the 
emitted energy and total angular momentum of the GWs 
into account; see, e.g., Eqs. (12) and (13) in [35]. Energy 
curves for SLy^-^^), SLy(^^\ and SLy*^°°^ are shown in 
Fig. 18 together with their pairwise differences (bottom 

panelb 


The differences between the spinning configurations 
and SLy(°°^ essentially quantify the repulsive spin-orbit 
(SO) interaction contribution to the binding energy, 
which is the dominant one for the dynamics. The energet¬ 
ics of SLy^-^^) are quite close to those of SLy^^^^. This 
happens because the leading order SO terms (oc L-Si/r^) 


The formulation of [ZQl uses non-spinning terms at the 3.5PN 
level, SO terms up to 4PN, spin-spin terms up to 2PN, and for 


the precession SO contributions up to the next-to-next-to-leading 
order are included. 
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t/M 

FIG. 17: norm of the Hamiltonian constraint (upper panel) 

and baryonic mass conservation (lower panel). Quantities are 
computed in level I — 1, i.e., the outermost Cartesian box of 
the numerical domain. 
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FIG. 18: Binding energy as a function of the reduced orbital 
angular momentum for the precessing simulation SLy^^^^ 
(blue dashed), SLy^°°* (solid orange), and SLy*-'’"''^ (solid 
green). The bottom panel reports differences between the 
curves. 


are identical in the initial data for and 

During the evolution, the difference between the SO inter¬ 
actions of SLy(^^) and SLy^-^^^ is solely due to the slight 
changes in the projection of the spins onto the angular 
momentum as they precess. For distances d ^ 16 — 8 this 
corresponds to corrections on the order of 10“'^ to the 
binding energy. The leading order spin-spin (SS) contri¬ 
bution {Ess = [3(n • 5'i)(n • S2) — {Si ■ 5'2)]/r'®, where 
n denotes the unit vector pointing from one star to the 
other) is exactly zero in the SLy*^^-^) initial data and or¬ 
der 10 “® in the SLy^^^^ initial data. During evolution the 
SS contribution of SLy^^"^) is of the order ~ 10 “®. This 
explains the differences between SLy^-^"^) and SLy^^^^ at 
the level of ~ lO^'^-lO”® due to these SO and SS correc¬ 
tions, which can be observed in Fig. | 18 [ However, sig¬ 
nificant differences between SLy*^^^) and SLy^^^^ can be 


The following discussion is based on leading-order PN expressions 
for the binding energy that can be found in, e.g., Eq. (2.7) of [106]. 


observed in some of the higher modes, notably the (2,1) 
mode, as we shall see in the following. 


2. Gravitational Waves 


We present the three largest amplitude modes of the 
GW signal in Fig. As in the nonprecessing case, the 
dominant emitter of GWs is the (2,2)-mode. However, 
interesting physical aspects are present in the subdomi¬ 
nant modes. The amplitude of the (2, l)-mode is modu¬ 
lated by the precession period, giving a second possibil¬ 
ity to extract the value of that period, which agrees with 
the estimate given above. Additionally, the amplitude 
of the modulation we observe is consistent with the ex¬ 
pected contribution to the (2, l)-mode from the binary’s 
dominant mass quadrupole radiation [which only appears 
in the (2, 2)-mode for nonprecessing systems] due to the 
precession of the orbital plane. In particular, since the 
initial angle between the total and orbital angular mo¬ 
menta is small in this case, Xl ■= E{L,j) ~ 0.05, de¬ 
creasing slightly over the evolution, as discussed in the 
previous section (see also Fig. 16), one can work to linear 
order in Al, where one finds that the maximum ampli¬ 
tude of the mass quadrupole’s contribution to the (2,1)- 
mode is 2Xl times the amplitude of the Newtonian mass 
quadrupole radiation.® One obtains the expression for 
the maximum mass quadrupole contribution to the (2,1) 
mode by noting that 2Ai is the largest angle the orbital 
angular momentum (originally aligned with the z-axis to 
a very good approximation) makes with the z-axis; cf. the 
PN expressions for the modes expanded in i in Eqs. (4.17) 
of |107j . (Note that t is defined as the angle between the 
orbital angular momentum and the initial total angular 
momentum, which they take to be along the z-axis in the 
coordinate system they use to define the mode decompo¬ 
sition.) Of course, there are contributions to the (2,1) 
mode from the current quadrupole, as well, but these are 
much smaller than the contribution due to precession, as 
is seen in Fig. since the current quadrupole contribu¬ 
tions are suppressed by a factor of vS compared to the 
Newtonian mass quadrupole radiation. 

The systems finally merges after ^ 30 GW cycles [in 
the (2,2) mode] at a frequency of ^ 0.128 (solid 

vertical line in upper left panel; for consistency with our 
previous work, we define the merger as the maximum of 
I/ 122 I) at t = 4:927M. 

Comparing with the other two simulations, 

we observe two main differences during the inspiral. First, 
without spin the merger happens earlier at t = 4839M, 
due to the fact that no repulsive spin-orbit interaction 
is present. In case of the SLy^^^^ simulation, the merger 
happens at t = 4960M, which agrees within At = 33M 


® Similarly, the amplitude of the (2, 2)-mode is the same as the 
amplitude of the Newtonian mass quadrupole radiation up to 
corrections that are suppressed by factors of vXj^S, or A|^, 
where u is the binary’s orbital velocity and (5 := — M^)/M ~ 

0 . 1 . 
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FIG. 19: The three most dominant multipoles of the GW {I, m) = (2, 2), (2,1), (3,3) for (left panels), (middle 

panels), and SLy*-®®^(right panels). We plot the real part of the modes and the dimensionless GW frequency in red. In the 
upper plots, the solid vertical lines mark the moment of merger (i.e., the maximum of |h 22 |), where a later merger in the two 
cases with spin is seen, due to the spin-orbit interaction. We observe that the (2,1) mode in the precessing case is dominated 
by the large contribution from the mass quadrupole at twice the orbital frequency due to the precession of the orbital plane, 
while it only has the much smaller contribution from the current quadrupole at the orbital frequency in the two nonprecessing 
cases. We extracted the wave at a radius of r = 446M. Some very small noise in the frequency at early times is visible and 
due to reflections at the boundary. After merger, the frequency calculation is less accurate and affected by larger oscillations, 
mostly unphysical. 


with the precessing simulation. The second observation 
is that the amplitude of the r/i 2 i is much smaller than 
for the precessing simulation. The non-zero amplitude is 
caused by the unequal masses of the two stars, but no 
clear imprint of the spin is visible (cf. the middle and 
right panels). Due to the small amplitude of the (2,1) 
mode in the nonprecessing simulations it is not very well 
resolved, which is clearly visible in unphysical sinusoidal 
oscillations of the frequency. The same holds to a lesser 
extent for the (3, 3)-mode. Additionally, the large fre¬ 
quency spikes present in the post-merger phase of the 
subdominant modes are caused by zeros of the amplitude. 

In order to further assess the differences between the 
SLy^-^^^ and SLy^^^^ waveforms, we compare the non¬ 
precessing data to the precessing ones after transforma¬ 
tion to the precessing frame. The transformation is per¬ 
formed with a method similar to the one used in binary 
black hole simulations in |53j . and the main result is 
shown in Fig.[^ In particular, the rotation of the '^4 im 
multipoles reads (Eq. (A9) in [55]) 

i 

^4im= , (5.3) 

m' = — l 


where d}^,^ are the Wigner d-matrices. We focus on the 
( 2 , 1 ) and ( 2 , 2 ) modes and perform the rotation as fol¬ 
lows: (i) we look for the Euler angles /3 ,7 for which | 4'4 21 1 
is minimal (Ref. [53] maximizes |^4 22 p + |'k 4 2 - 2 P); (h) 
we split the dataset in several chunks (vertical dashed 
lines in Fig. [20] ) and fit the obtained Euler angles with 
low-order polynomials. This step optimizes the fits and 
minimizes numerical oscillations. We do not need to 
worry about the third Euler angle a here, as it is irrele¬ 
vant in the case we consider, where we only look at the 
magnitude of the modes of ^ 4 . From Fig. W one sees that 
(i-®-; version in the nonprecessing frame) is 


almost equivalent to 


(n) 


The effect is very clear 


in the (2,1) mode. This preliminary result suggests it 
will be possible to model precessing BNS waveform using 
aligned spin BNS models for moderate spin magnitudes 

(cf. [55]). 


The merger remnant is a hypermassive neutron star 
(HMNS), which mostly emits in the (2, 2) channel at early 
times. After t « 5900M the amplitude of the (2, 2) mode 
decreases until the (2,1) and (3, 3) modes have the same 
amplitude as the (2,2) mode. Note that at this time 
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FIG. 20: Amplitude of the (2,2) and (2,1) modes of \l /4 for 
the and simulations, along with the same 

amplitudes for SLy^^^* in the nonprecessing frame (denoted 
by \[' 4 ) obtained as in [53]. Note that we only consider the 
inspiral here. The transformation is done in chunks (separated 
by vertical dashed lines). The agreement between 
and I after transformation to the non-precessing frame 

is clearly visible. The noise around u « 1500M is caused by 
reflections of the outer boundary. 


the (2,0) and (4,4) modes also have comparable am¬ 
plitudes. We find a frequency shift of the / 2 -frequency 
[the dominant frequency in the (2, 2) mode] of ~ 60 Hz 
due to the additional angular momentum of the HMNS 
formed by the spinning configurations; the origin of such 
frequency shifts was discussed in detail in |3S]. The es¬ 
timated / 2 -frequencies are 2.75, 2.79, and 2.81 kHz, for 
SLy(°°), SLy(^^), and SLy^-^^^ respectively. 


C. Effect of eccentricity reduction on waveform 
phasing 

Although an eccentricity reduction procedure for BNS 
quasicircular initial data has been already presented 
in [29j . its performance on the relevant observable quan¬ 
tity, i.e., the GW phase and amplitude, has not been eval¬ 
uated directly. Here, we investigate the effect of eccen¬ 
tricity reduced data on the GW phasing and amplitude 
by a direct comparison with data in which eccentricity 
reduction has not been performed. Such a comparison 
is particularly important since the eccentricity reduction 
procedure is computationally expensive and might be not 
necessary for certain applications or when the data are af¬ 
fected by larger uncertainties due, for example, to trun¬ 
cation errors. 

We compare the GW phase of two SLy q = I runs: 
In one case we evolve initial data with Bd = 1.241 x 10“^ 
(Iter 0) and in the other data with = 8.7x 10“"^ (Iter 3); 
see Table UTl 

We focus on the £ = m = 2 multipole and omit the sub¬ 
script in the waveform quantities. Waveforms are aligned 
on the interval [ti,t 2 ] = [1000,6000] ~ [370M, 2222M] 
shifting by constant time and phase offsets T, 4). The 



FIG. 21: GW phasing of the evolved initial data without 
(Iter 0) and with (Iter 3) eccentricity reduction for the SLy 
EOS ? = 1 configuration. Vertical dotted lines denote the 
region over which we align the waveforms. 

latter are determined by minimizing the function |108j 
G(r, $) = / ' |<^i(t) - Mt + T)- (5.4) 

Jti 

where ^i _2 denotes the GW phase of the two datasets. A 
more robust alignment procedure based on a frequency 
interval can also be used but the current procedure 

is sufficient for our purposes. 

We present the results in Fig. The phase difference 
A 022 oscillates between [—0.06,0.06] rad during the ^ 21 
GW cycles, and it is essentially flat up to merger, t ~ 
2780M. Similarly, the amplitude of the non-eccentricity 
reduced data (Iter 0) oscillates around the eccentricity 
reduced ones (Iter 3); the amplitude oscillations are ~ 
5% at early times t ~ 370M and decrease as the system 
approaches merger. 

Overall, these results show that the use of eccentricity 
reduced data with Bd ~ 10“^ (about 3 iterations of our 
procedure) improves the waveform quality for GW model¬ 
ing purposes, and should be employed in future precision 
studies of the gravitational waveform. However, eccen¬ 
tricity reduction is likely to be effective only if combined 
together with an improvement of other source of errors, 
notably truncation errors. We notice in this respect that 
A(j> ~ 0.12 rad is at least a factor two smaller than the 
typical uncertainty introduced by truncation errors at the 
resolutions employed here |100j . 

VI. CONCLUSIONS 

Due to advances in the construction of constraint 
solved and consistent initial data, simulations of binary 
neutron stars in full general relativity are now able to 
cover more of the binary neutron star parameter space ac¬ 
curately. In particular, it is now possible to study differ¬ 
ent EOS large mass ratios [Ml 177] . spinning neu¬ 

tron star configurations [5H], highly eccentric setups HZ], 
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and neutron stars on orbits with less residual eccentricity 
than the standard ones 1109] . We have recently up¬ 

graded the SGRID code to be able to generate consistent, 
constraint-solved initial data with the ability to vary all 
these parameters.The most noteworthy improvement 
was the combination of the constant rotational velocity 
approach (from [23) 1 with a generic specification of the 
symmetry vector (from |17j l allowing arbitrary eccentric¬ 
ities (including eccentricity reduction). 

In this paper, we have exhibited SGRID’s ability 
to generate binary neutron star initial data for many 
situations of interest and included dynamical simulations 
of some generic configurations evolved with the BAM 
code. 


A. Quasi-equilibrium sequences 

We have constructed the following quasi-equilibrium 
setups: 

Binaries in the constant rotational velocity approach. 
With the help of quasi-equilibrium sequences we studied 
spinning neutron star configurations for different equa¬ 
tions of state and characterized the spin-orbit contribu¬ 
tion to the binding energy. Our results show that the 
spin-orbit interaction can be well approximated by post- 
Newtonian theory within the uncertainty of our numerical 
method. 

Highly eccentric binaries. We also constructed highly 
eccentric sequences following the description of but 
with the advantage of solving the elliptic equation for 
the velocity potential, which results in smaller artificial 
density oscillations by a factor of ^ 5. This will allow 
a more detailed analysis of tidally-induced oscillations in 
the neutron stars than was performed in |30j using su¬ 
perposed initial data, which led to much larger initial 
oscillations. Additionally, we compared our eccentric se¬ 
quences with PN calculations in various ways, where our 
analysis showed that, as expected, close agreement is only 
obtained in the limit of large separations or small eccen¬ 
tricities. 

Eccentricity reduced binaries. We can use the same 
technology that allows us to create highly eccentric orbits 
to reduce the eccentricity present in standard binary neu¬ 
tron star initial data constructed using a helical Killing 
vector. Here we iteratively adjust our eccentricity and 
radial velocity parameters, similar to the procedure pre¬ 
sented in P5] . 

Varying compactnesses and mass-ratios. As additional 
parameters, we considered different mass ratios and 
compactnesses. We were able to construct a sequence for 
a mass ratio oi q = 2.06 and equal-mass binary neutron 


Note that here we only consider the “inspiral parameter space,” 
consisting of the purely relativistic hydrodynamic parameters (ec¬ 
centricity, masses, spins, and EOS) that can be measured with 
an inspiral gravitational wave signal with current or proposed 
gravitational wave detectors in physically expected scenarios. 


stars with compactness up to C = 0.23. 

This shows that with its recent upgrades, SGRID al¬ 
lows one to construct binary neutron star configurations 
in a substantial portion of the possible inspiral parameter 
space (and one can vary all the relevant parameters in¬ 
dependently). Of course, there is a considerable amount 
of physics that we do not include here. However, the 
physics we have neglected does not affect the inspiral at a 
level that can be detected via gravitational waves for the 
parameter values expected to be present in neutron star 
binaries. Nevertheless, this missing physics can still play 
an important role in the merger or potentially produce 
other interesting effects, and includes magnetic fields, 
elasticity (e.g., in the solid crust), and composition, all 
of which would need to be appropriately incorporated 
in the initial data. Of this additional physics, the 
inclusion of magnetic fields is likely the most pressing. 
Unfortunately, there is no known method for including 
magnetic fields consistently in constraint-solved initial 
data: All simulations of magnetized binary neutron 
stars (e.g., HUT]) add the magnetic field by hand after 
constraint solving. Developing such a method would 
be a useful advance in binary neutron star initial data 
construction. 


B. Dynamical Evolutions 

To ensure that the constructed data are suitable for 
dynamical simulations, we have evolved three configura¬ 
tions. 

q = 2.06 run. As a first example, we evolved the high¬ 
est mass ratio ever considered in a full general relativis¬ 
tic binary neutron star configuration. The configuration 
consisted of a, q = 2.06 setup with the MSlb EOS. Be¬ 
cause of the high mass ratio and the rather stiff EOS, 
we observed mass transfer between the two stars several 
revolutions before merger. During this process material 
with a rest-mass of ^ (2 — 3) x 10~^Mq is transferred, 
with an average accretion power of ^ 10®^ erg s“^. Dur¬ 
ing the merger process, ~ 7.6 x get unbound 

and are released from the system with a kinetic energy 
of ~ 4 X 10®° erg. The ejecta process happens primarily 
due to torque in the tidal tail of the lower massive star 
and forms a spiral like pattern. Due to this anisotropic 
ejection of material, the merger remnant receives a large 
kick of 0(100) km s~^. The final merger remnant can 
be characterized as a supramassive neutron star, which 
is not expected to collapse on dynamical timescales. An 
investigation of the merger remnant’s GW spectrum, in¬ 
cluding more than just the dominant (2, 2)-mode, reveals 
that many of the peak frequencies are harmonically re¬ 
lated to high accuracy. 

Processing and unequal masses run. As a second ex¬ 
ample, we considered the first precessing binary neutron 
star merger simulation. Contrary to most BNS investiga¬ 
tions (except, e.g., [6^ Ifl0| )> we present more than just 
the dominant (2, 2) mode and find a clear imprint of the 
precession in the subdominant (2,1) mode, where the am- 
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plitude is modulated by the precession frequency. Consid¬ 
ering for comparison a simulation with the same leading- 
order spin-orbit interaction, we show that the relation 
between the gauge-invariant binding energy vs. reduced 
orbital angular momentum up to the merger exhibits only 
a minor imprint of the precession for the spin magnitudes 
we consider. Regarding the post-merger GW spectrum, 
we observed a frequency shift of the dominant /2 mode 
due to the spins of the binary’s components (cf. [55]1. 

Reduced eccentricity run. Finally, we have performed 
simulations of an equal-mass configuration with and with¬ 
out eccentricity reduced initial data and we have quan¬ 
tified the differences in the waveform’s amplitude and 
phase. We found that although the eccentricity reduc¬ 
tion improves the waveform quality, one also needs to 
reduce other errors in the waveforms, notably truncation 
errors, in order for the improvement due to eccentricity 
reduction to be effective. 
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Appendix A: Astrophysical predictions for more 
extreme masses, mass ratios, and spins for binary 
neutron stars 

Here we assess the prospects for binary neutron star 
mergers with more extreme masses and spins—such as 
those we simulated in Sec. |V] —actually occurring in na¬ 
ture. (Note that this assessment only covers the masses 
and spin magnitudes, not the equation of state or spin 
misalignment.) 


1. Large and small masses and larger mass ratios 

There is good evidence that the population of neutron 
stars in the universe extends at least from ~ IAIq up 
to 2Mq. In particular, there are two precise measure¬ 
ments of high-mass (~ 2Mq) neutron stars in binaries 
with white dwarfs [63l |64] as well as some less precise 
measurements of low-mass {IMq) neutron stars in X-ray 
binaries mm- There are also recent measurements of 
low compactnesses for isolated neutron stars m, which 
would imply quite small masses ^ IAIq for many of the 
equations of state we consider.Additionally, there is 
recent work that suggests that the initial (pre-accretion) 
mass of the low-mass millisecond pulsar J0751-I-1807 
(whose present mass is 1.26 ± O.I 2 M 0 ) could have been 
as low as I.IMq m- 

The theoretical bounds on the neutron star mass al¬ 
low for an even larger mass range (and thus, in princi¬ 
ple, large mass ratios, up to ^ 3), as some EOSs have 
maximum masses of ^ 3M0 (e.g., the ~ 2.8Mq maxi¬ 
mum masses for the MSI and MSlb EOSs we consider 
in this work; see Table |T|. However, the minimum mass 
of neutron stars in the universe is likely around the mini¬ 
mum observed mass of ^ IMq. While the minimum mass 
of a star constructed from a cold dense matter equation 
of state is quite small (< O.IM 0 ), the minimum mass 
of a hot protoneutron star is considerably larger, 0.89- 
I.I 3 M 0 for the models considered in |116) . This mini¬ 
mum mass provides a practical lower bound on neutron 
star masses formed from supernovae, barring formation of 
lower-mass stars by fragmentation (see, e.g., (nil), which 
is quite speculative. Moreover, the high kicks that neu¬ 
tron stars formed by fragmentation would be expected to 
receive make them unlikely components of binaries. Ad¬ 
ditionally, there is a further restriction from the baryonic 
mass of the iron core of the supernova progenitor, which 
gives a minimum mass of 1.15-1.2M0, as discussed in 
Sec. 3.3 of [T], though there are uncertainties in both 
these bounds due to uncertainties in supernova physics. 
(Note that Tauris, Langer, and Podsiadlowski |118) esti¬ 
mate that the minimum mass of a neutron star formed 
in an ultra-stripped supernova is I.IM 0 .) See [HU] for a 
general review of neutron star masses. 


Note that these authors quote compactnesses in units of Mq /km, 
not the dimensionless compactnesses we use in this work, as is 
mentioned explicitly in the caption to Fig. 1 of nni- 
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As mentioned in Sec. |Tj the mass range of the observed 
binary neutron star systems is much smaller, particu¬ 
larly if one only considers the six systems that will merge 
within a Hubble time, where the minimum and maxi¬ 
mum masses are 1.25 and 1.44Mq (PSRs J0737—3039B 
and B1913-I-16, the less massive star in the Double Pul¬ 
sar and the Hulse-Taylor pulsar, respectively), and the 
largest observed mass ratio is 1.07, for the Double Pul¬ 
sar (see, e.g., Table 1 in pQ [5] and Table 3 in |119j l. 
Thus, one might naively not expect to see a very wide 
range of masses (and thus mass ratios) in many binary 
neutron star coalescences. However, for certain values of 
poorly constrained parameters, population synthesis cal¬ 
culations (e.g., the Synthetic Universe models from Do- 
minik et al. m) predict the existence of binary neutron 
stars (formed “in situ," i.e., not by dynamical capture) 
with masses over the entire observed range of neutron star 
masses, and some systems with reasonably large mass ra¬ 
tios. 

Let us now consider the predictions of the Synthetic 
Universe population synthesis data available online [120) ; 
these are the standard model and Variations 1-15 of Do- 
minik et al. [7] , where each of the variations varies one of 
the poorly constrained parameters in the calculation—see 
Table 1 in [7] . (See |121] for a study of the predicted grav¬ 
itational wave detection rates using a few of these models 
and m for a study of the effects of varying certain ini¬ 
tial conditions.) Additionally, each of these models has 
four further variants, given by the four combinations of 
two choices for the metallicity (solar and 0.1 solar)and 
two treatments of the common envelope phase of the bi¬ 
nary’s evolution (submodels A and B, which correspond 
to the optimistic and pessimistic predictions, respectively, 
for the fate of binaries which enter the common envelope 
phase when the donor is in the Hertzsprung gap). These 
models all assume a minimum neutron star mass of IMq 
(as mentioned in m) and a maximum mass of 2.5Mq 
(except for Variations 5 and 6, which assume a maximum 
mass of 3 and 2Mq, respectively). 

All of these models predict a galactic binary neutron 
star merger rate above the estimated lower bound of 
2.1 Myr“^ inferred from observations (see the discussion 
in Sec 4.1 of |124| 1 at solar metallicity, except for Varia¬ 
tion 1 (both submodels), and submodel B in Variations 2, 
4, and 12 (see Table 2 in H). We still show results 
for these models, for comparison, particularly since that 
lower bound is not particularly firm. 

We show the total number of coalescing binaries in each 
of these models in Fig. marking the numbers of sys¬ 
tems with individual masses > 1.5, 1.75, and 2Mq. We 
also show the maximum and minimum individual masses 
and mass ratio present in the coalescing systems in that 
figure. (Here we select the systems that coalesce within 
10 Gyr from the formation of the binary, the criterion 


Note that the metallicity of objects in the universe varies, gener¬ 
ally increasing for more recent formation times, as discussed in, 
e.g., [na. These two choices of metallicities are intended to give 
an indication of the effects of metallicity on these calculations. 





a A02 a A002 a B02 a B002 


FIG. 22: The total number of coalescing binary neutron star 
systems in the various Synthetic Universe population synthe¬ 
sis models, with the numbers with individual masses greater 
than > 1.5, 1.75, and 2Mq also marked (top), along with the 
maximum and minimum individual masses (middle) and max¬ 
imum mass ratio (bottom) present in these models. Here “S” 
denotes the standard model and “Vn” denotes the nth varia¬ 
tion. The four variants are denoted by different colored bars, 
where the letter (A or B) gives the submodel and “02” and 
“002” correspond to solar metallicity and 0.1 solar metallic¬ 
ity, respectively, the notation used on the Synthetic Universe 
webpage. In the top plot, the number of coalescing systems 
with individual masses > 1.5, 1.75, and 2Mq are marked with 
dense hatching, less dense hatching, and dotted hatching, re¬ 
spectively. 
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FIG. 23: A histogram of mass ratio of coalescing systems for 
the various Synthetic Universe models that have a maximum 
(coalescing) mass ratio of > 1.7. The number of coalescing 
systems in each model is given after the model’s name (with 
the same notation as in Fig. 221. Each bin has size 0.1 and is 
labeled by the maximum mass ratio present in the bin. The 
models are ordered to help clarify the plot. 


for the population of potentially observable coalescing 
binaries used in [7].) We see that many of the models 
predict at least some coalescing systems with individ¬ 
ual masses close to 2Mq (at least above 1.75Mq). The 
minimum individual masses present are mostly close to 
I.IM 0 , though in a few models they approach the mini¬ 
mum mass of IMq, while in Variation 10 (which assumes 
a delayed supernova engine), the minimum mass in the 
low metallicity case is 1.2Mq. The maximum total mass 
is never more than 20% smaller than twice the maximum 
individual mass, and in almost half of the cases is within 
5% of it. 

The minimum total mass is always very close to 2.2Mq, 
except for Variation 1, where it is 2A9Mq (for all four 
cases), and Variation 10, where it is between 2.3 and 
2.5Mq, depending on the case. Indeed, mergers with 
a stable neutron star as the final remnant may even be 
rather common, if the neutron star maximum mass is 
> 2.5Mq as is assumed in all but Variation 6: One finds 
hundreds to thousands of coalescing systems with total 
masses below this limit in all but a handful of these Syn¬ 
thetic Universe models. Moreover, these numbers are a 
significant fraction of the total number of systems in each 
of the population models (in a number of cases well above 
90%, including the standard case with solar metallicity 
and submodel B, where the fraction is 99.88%). There 
are even tens to hundreds of coalescing systems with to¬ 
tal masses less than 2.25Mq in many models. In the 
models from m, which vary various initial conditions 
using a base model quite similar to the standard model 
from [ 7 ], one does not find many of the more extreme 
systems we consider here, since these models do not in¬ 
clude the choices for the binary evolution parameters that 


generate such systems. 

Since there is a wide range of individual masses in these 
models, one might expect there to also be a wide range 
of mass ratios present, which we indeed find to be the 
case, as illustrated in the bottom panel of Fig. We 
also show a histogram of the number of systems with 
different mass ratios for the models that have maximum 
mass ratios > 1.7 in Fig. |23| The maximum mass ratio 
in the full (not just the coalescing) population is 1.94, 
found in Variation 10 with solar metallicity and either 
submodel. The maximum mass ratio in the coalescing 
population is 1.86 and is also found in Variation 10 with 
solar metallicity, though here just for submodel A. Note 
that this variation assumes a delayed supernova engine, 
and thus may be unphysical, as discussed in [7] . We there¬ 
fore note that Variation 7 (which assumes low supernova 
kicks) with solar metallicity gives a maximum mass ra¬ 
tio of 1.83 in both the full and coalescing populations for 
both submodels, as does Variation 14 (which assumes a 
weakly bound common envelope) for 0.1 solar metallicity 
and submodel A. 

While the very largest mass ratios are indeed uncom¬ 
mon even in these more extreme population models (see 
Fig. 23), two models predict > 50 coalescing systems 
with mass ratios > 1.7 (Variation 4, which assumes a 
very weakly bound common envelope, with submodel A 
and either metallicity). These models are on the extreme 
side, but even the standard model predicts a single mass 
ratio 1.5 coalescing system (with submodel A and solar 
metallicity). 


2. Larger spins 

Compared to the spins of millisecond pulsars (where 
the largest spin known is 716 Hz |125j ). the members 
of known binary neutron stars are not spinning nearly 
so rapidly. The shortest spin period observed so far 
is 22.7 ms (corresponding to a frequency of 44 Hz) 
for the more massive star in the Double Pulsar, PSR 
J0737—3039A, which is likely not to spin down too much 
by the time the system merges, as discussed in Ap¬ 
pendix |A^ PSR J0737—3039A has a dimensionless spin 
of j G [0.02, 0.03], where the uncertainty comes from the 
uncertainty in the EOS. 

However, for a direct estimate of the typical spin in bi¬ 
nary neutron star systems, we are restricted by the fact 
that we currently only know a small sample of the bi¬ 
nary neutron star systems in our galaxy (around twelve). 
Furthermore, the existence of the sizeable population of 
rapidly spinning neutron stars, with 207 known pulsars 


While PSR J1807—2500B (NGC 6445B) has a spin of 239 Hz, it 
is at present unclear whether its companion is a neutron star or 
a white dwarf (The uncertainty about the nature of the 

companion is somewhat larger for this system than it is for most 
binary pulsars. This is also the most massive companion of a 
fully recycled neutron star.) Additionally, this system will not 
merge within a Hubble time, so it would not contribute directly 
to binary neutron star merger rate calculations. 
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with spins over 200 Hz |127] suggests that there should 
be a population of binary neutron stars where at least one 
star has a significant spin at merger. (Note that 200 Hz 
corresponds to a dimensionless spin of ~ 0.1 for the EOSs 
considered in this work.) This would most likely be the 
heavier star, which could have its spin increased by ac¬ 
cretion from its companion when its companion is still 
a post-main sequence star, a process known as recycling 
(see, e.g., |129) '). (The recycling process reduces the star’s 
external magnetic field, which will likely allow the rapid 
spin obtained from accretion to persist until merger.) 

However, one needs to accrete a fair amount of mat¬ 
ter in order to spin a neutron star up to high frequencies 
(^ O.IMq for frequencies of ^ 500 Hz is quoted in Sec. 7 
in m\), and one does not expect to accrete lots of matter 
when the companion is a neutron star progenitor. In par¬ 
ticular, Tauris, Langer, and Podsiadlowski |118j consider 
the ultra-stripped supernova binary neutron star forma¬ 
tion channel, which they claim is the primary channel 
for forming binary neutron stars that will merge within 
a Hubble time. In their calculations, they obtain a maxi¬ 
mum spin of ~ 40 Hz, assuming Eddington-limited accre¬ 
tion onto the neutron star, and ~ 90 Hz, assuming three 
times Eddington accretion (which they mention is eas¬ 
ily obtained, and for which there may be evidence in the 
inferred formation channels of other systems); see their 
Sec. 6.1.1. 

Additionally, MacLeod and Ramirez-Ruiz |13()j find 
that even though the neutron star only accretes < O.IMq 
during the common envelope phase in their simulations, 
the neutron star could still be spun up to ~ 250 Hz in 
the cases where one accretes close to O.IMq. Of course, 
if the neutron star experiences a significant enough spin- 
down after the accretion episode, then these high spins 
will reduce substantially before merger. However, one 
expects the magnetic field to be reduced by accretion, as 
mentioned above, and for reasonably small periods and 
magnetic fields after the common envelope phase (e.g., 
^ 4 hours and ^ 10® G), the spin at merger can still be 
^ 200 Hz. See the next subsection for further discussion 
of the issue of determining the spin at merger. 

One can also imagine forming a double neutron star 
with a member with a high spin through dynamical for¬ 
mation via binary-single (or even binary-binary) interac¬ 
tions in dense stellar regions, such as globular clusters 
(see, e.g., 0 HH]), where one swaps out the companion 
that recycled the highly spinning neutron star. This pro¬ 
cess is discussed as a likely formation channel for the 
binary containing the 239 Hz pulsar J1807—2500B |126) . 
which resides in a globular cluster, and whose compan¬ 
ion may be another neutron star. There is even the ex¬ 
otic possibility of forming a binary where both neutron 
stars have millisecond periods through this channel, or 
through the double recycling scenario proposed by Sig- 
urdsson and Hernquist [9], where the neutron stars’ main 


This is out of 2525 known pulsars. Also note that 20 out of the 
23 known pulsars in the rich globular cluster 47 Tucanae have 
spins over 200 Hz 


sequence companions are disrupted and this material re¬ 
cycles both neutron stars. 


3. Neutron star spin predictions at merger 


In addition to considering the purely theoretical 
prospects for relatively high spin in merging binary neu¬ 
tron stars, one can also consider the spins at merger of 
known pulsars whose companion is a neutron star. The 
simplest way to estimate the spin at merger is to assume 
that the pulsar’s observed spindown is due solely to mag¬ 
netic dipole radiation. This was already done for the 
fastest-spinning neutron star known in a neutron star bi¬ 
nary (J0737—3039A, the more massive star in the Dou¬ 
ble Pulsar) by one of us in [53]. There it was found that 
this pulsar’s spin will only decrease slightly, from 44 to 
37 Hz, in the 85 Myr from now until the system merges, 
with this assumption. Although the observed spindown 
is likely not due entirely to magnetic dipole radiation, we 
will show that the estimate for this particular pulsar is 
valid without this assumption (with reasonable assump¬ 
tions about the size of the pulsar’s braking index). 

For magnetic dipole radiation the pulsar’s braking in¬ 
dex (n := where v is the pulsar’s spin frequency) 

is 3. However, all reliably measured values of the braking 
index (only 8) are less than 3 (and can be as low as ~ 1); 
see, e.g.. Table I in m- Note that all of the pulsars 
with reliably measured braking indices are much younger 
than J0737—3039A, with ages of at most ~ 10"^ years, 
as opposed to J0737—3039A’s age of ~ 10® years. If one 
looks at these braking indices versus the pulsars’ char¬ 
acteristic ages from the ATNF catalogue , one sees 
that the significantly older pulsars all have lower—but 
more uncertain—braking indices than the younger ones. 
(See Fig. 7.5 in |132j for an illustration, albeit without 
error bars. They also plot less well-determined braking 
indices for much older pulsars, finding quite large values 
up to ~ 10®, though they are also likely quite uncer¬ 
tain.) All these pulsars are also more slowly rotating 
than J0737—3039A, with a largest frequency of ~ 30 Hz 
for the Crab pulsar. Nevertheless, it turns out that the 
prediction of J0737—3039A’s spin at merger is quite in¬ 
sensitive to the braking index assumed, since the spin- 
down timescale is considerably longer than the time to 
merger. 

Specifically, the period evolution of a pulsar with a 
generic (constant) braking index n is given by P = 
where AT is a constant. We can solve this in 
terms of the pulsar’s period at t = 0, Pq, and express 
K in terms of the pulsar’s characteristic age at t = 0, 
To ■■= Po/(2Po), giving 


P{t) = Po 


= Po 



n — 1 t 


l/(n-l) 


2 ro_ 

It 2 — n l't \ 

2 To 8 \To) 


2 


+ 0 



(Al) 


Here we have expanded to second order in t/ro to illus- 



















32 


trate that the correction term due to the braking index is 
quite small when t/rg is small. In case of J0737—3039A, 
where tq = 210 Myr, so tmerger/T'o — 0.4, the corrections 
due to possible deviations of the braking index from 3 
are < 5% (considering n G [0,5], where n = 5 corre¬ 
sponds to pure gravitational wave damping). Addition¬ 
ally, for the very large magnitude (and likely quite uncer¬ 
tain) braking indices found for pulsars of about the age of 
J0737—3039A discussed above, the change in the pulsar’s 
spin until merger is negligible. 

Of course, this calculation still assumes a constant 
braking index and K. If, for instance, J0737—3039A has 
a large buried magnetic field from the accretion episode 
that is thought to have spun it up (its external field is 
only 6.3 x 10® G, while its companion’s is ^ 10^® G) |133j . 
and enough of this field becomes unburied before the bi¬ 
nary coalesces, then this could spin the pulsar down far 
more than is predicted by the calculation above: Since 
the spin-down rate goes as the square of the magnetic 
field, an increase in the external field by an order of mag¬ 
nitude would decrease the characteristic age r by a factor 
of 100, so one could easily obtain a final spin that was 
smaller by a factor of 2 or more if this larger external 
field is present for ^ 10 Myr or longer. And if a ~ 10^^ G 
field emerges completely, r will decrease to ~ 10^ years, 
so such an external field could spin the star down by an 
order of magnitude in only ^ 1 Myr. 

However, this scenario relies on the pulsar’s having ac¬ 
creted a small enough amount of matter to allow a signifi¬ 
cant amount of field to become unburied in the remaining 
time before merger. As illustrated in, e.g.. Fig. 1 of |134j . 
the unburial timescale is very sensitive to the amount of 
accreted material (with factors of 2 leading to order of 
magnitude changes). Since there is no way of which we 
are aware of estimating the amount of material accreted 
onto J0737—3039A to within a factor of 2, and the en¬ 
tire picture is further complicated by the possibility of 
magnetic field decay (e.g., Ohmic decay), not just burial 
(see, e.g.. Sec. 10.3.1 in |135j for a brief review), we do 
not pursue this possibility further here. 


Appendix B: Alternative derivation of a first integral 
to the Euler equation 


where in the last step we have used Eqs. ( 2.19) an d (2.201. 

Using u ■ u = —1, the Euler equatit m (|2.11 ) can be 
written in the Carter-Lichnerowicz [13611137j form: 


u ■ dp = 0, 


Using Eq. (2.19), Eq. (B4) can be rewritten as 

k ■ dp -b V ■ dp = 0. 


(B4) 


(B5) 


We now rewrite k ■ dp using the Cartan identity (B1) to 
obtain 


CkP - d(fc • p) -(- y • dp = 0, 


(B6) 


where CkP vanishes by assumption if the symmetry vec¬ 
tor k Lie-derives the flow. We could then obtain a first 
integral k ■ p = const for irrotational (dp = 0) or coro- 
tational {V = 0) flow, since then the last term is zero as 
well. In the general case of spinning neutron stars, how¬ 
ever, that is not true. Nevertheless, we can again make 
use of our assumptions ( 2.24b[ ), (2.24c), and (2.27) (this 
time without projecting onto the three dimensional slice), 
i.e., 


£fed()) « « 0, (B7) 

where the latter two can be combined to get 

£ p u; = -b ic « 0. (B8) 

pO pB 


Now it is possible to rewrite the first and the last term 
of equation (B6) with the goal of finding a first integral. 
Therefore, we can write C^p — £fc(d^ -b to) « to 
simplify the first term. For the last term, we will apply 
Cartan’s identity once again and the definitions from (B2) 
to write 


V ■ dp = y • dm = Cvw — d(y • w). 


(B9) 


Substituting these two terms into (B6) and exploiting 
linearity of the Lie derivative yields 

d{k p + V ■ w) Ck+vw = £^_|_ p, m « 0. (BIO) 

p^ 

which gives rise to an approximate first integral 

k ■ p + V • m « const (Bll) 


We have already outlined a possible derivation of a first 
integral to E uler’s equation for rotating neutron stars, 
see Sec. BA and [23l|24j. In this appendix we obtain the 


first integral using a much shorter derivation based on the 
Gartan identity, which relates the Lie derivative operator 
£ to the exterior derivative operator d. In particular, for 
a differential form a; and a vector u one has 


CuUj = u ■ do) -b d(M • u)). 


(Bl) 


where a dot denotes contraction between adjacent indices. 
We will now use the canonical momentum 


p = hu = d(j) ■ 


(B2) 


as introduced in Eqs. (|2.20 |) and (|2.23|), together with the 

(B3) 


additional definition (2.25), to write 

p = -b m = p°(fe -b y), 


In a final step, using the normalization condition u ■ u = 

— I, one can straightforwardly show that fc • p -b y ■ w = 

— ^ — y • d(j), which corresponds to equati on (|2.30 ). 
Notice, however, that the assumptions ( |B7[ ) used in 

the derivation here are slightly stronger than the origi¬ 
nal assumptions (2.24b), (2.24c) and (2.27), which make 


assumptions only about projections onto the three dimen¬ 
sional slice. 


Appendix C: Single CRV-stars 

1. Comparison with rigidly rotating stars 

As explained in [53] , the pa rticul ar choice of the angu¬ 
lar velocity given by Eq. (2.35) leads to a negligible 
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shear, thus any substantial differential rotation can be 
neglected. We provide corroborating evidence for this 
statement here. For this purpose we construct single ro¬ 
tating neutron stars in the CRV approach and compare 
them with rigidly rotating stars. We compute the rigidly 
rotating stars with the project Nrotstar of the publicly 
available LORENE library m- Nrotstar solves the Ein¬ 
stein equations with a self-consistent field method and 
multi-domain spectral methods. To compute CRV stars, 
we use SGRID and the formalism described in Sec. m 
In particular, we take the approximate symmetry vec¬ 
tor k to be the timelike Killing vector dt- We also use 
UA = ns = 26, = 8 , and ncart = 22 . 

In Fig. 24 we present sequences of CRV and rigidly ro¬ 
tating stars with the simple polytropic r2.72 EOS for dif¬ 
ferent central enthalpies (i.e. different baryonic masses). 
There is no evidence that the result would be different 
for piecewise polytropes. 

Due to the fact that in the CRV approach does not 
agree with the frequency measured by an observer at in¬ 
finity (as already outlined in [55])) we compute CRV con¬ 
figurations for a constant and find a corresponding se¬ 
quence of rigidly rotating stars by choosing the frequency 
in Nrotstar such that both methods give similar re¬ 
sults. We obtain the SGRID data by varying the baryonic 
masses Mb G [1.1,2.0] with a spacing of AMb = 0.1. We 
clearly see that for frequencies below 300Hz, rigidly ro¬ 
tating data and CRV data agree well. However for larger 
angular momentum larger discrepancies occur. This can 
be caused by (i) th e reduced accuracy for faster rotating 
stars (see Sec. |C 3| ) and (ii) the fact that the members of a 
sequence with different masses and the same = const 
may not all correspond to the same observable frequency. 


2. Empirical a;-j-Relation 


In addition to the comparison with rigid rotating stars, 
we want to answer the question of how the dimensionless 
spin of a single CRV-star can be obtained from input 
parameters of the SGRID code, namely the EOS, the 
baryonic mass Mb, and the angular velocity vector w®. 
A phenomenological model to obtain a rough estimate 
of the star’s spin value would reduce the computational 
costs to find accurate initial data for particular configu¬ 
rations. Eor this reason we choose four EOS: SLy, ALF2, 
MSlb, and r2 with baryonic masses in the range [I.l, 1.7], 
spanning a range in the compactness of C G (0.09, 0.20). 

The dimensionless angular momentum j of a neutron 
star is given by 



^central 


FIG. 24: Comparison of single neutron stars with rigid rota¬ 
tion and CRV rotation using a simple polytropic EOS (r2.72). 
The values for in the CRV-approach are 0.0, 0.005, 0.01, 
and 0.015. (These data are shown as crosses.) The solid lines 
are computed for rigid rotating stars with the LORENE li¬ 
brary. 


in the form 


j = f{C,Mk)uj. (C2) 

Also, the gravitational mass of the single star (i.e., 
Madm) for this spacetime is not known in advance and is 
thus absorbed in the function /. We find with numerical 
experiments the following expression: 

jut — oi (1-|-wiMb)(1-l-ciC-fC 2 C^-|-C 4 C^)tu, (C3) 

where the parameters oi = 88.8131, toi = 1.39522, 
Cl = -19.003, C 2 = 152.99, C 3 = -570.678, C 4 = 
806.896 are obtained from a fit of the data presented 
in Eig. and where we use the compactness C of an 
irrotational star with the same baryonic mass for sim¬ 
plicity. Specifically, we computed five different values 
oj = (0.000,0.002,0.004,0.006,0.008) for each of the 
baryonic masses (1.1,1.2,1.3,1.4,1.5,1.6,1.7) and for the 
EOSs SLy, ALE2, MSlb, and F2. We want to stress that 
the resulting relation is just empirical and probably does 
not represent any underlying physical properties. 


JaDM _ li^ohs 
~ m2 ’ 


(Cl) 


3. Rapidly rotating neutron stars 


where Jadm and Madm are the spacetime’s ADM angu¬ 
lar momentum and mass, respectively, and I the moment 
of inertia of the star. Now Wobs) the rotational period an 
observer at infinity would measure, is not known a priori, 
but probably depends linearly on the angular velocit y a; , 
for slowly rotating neutron stars. We thus recast Eq. (Cl I 


The fastest spinning neutron star observed so far is 
PSR J1748—2446ad, with a spin period of 1.4 ms, cor¬ 
responding to a frequency of 716 Hz |125j . This corre¬ 
sponds to a dimensionless spin of j G [0.3, 0.6] (for the 
EOSs given in Tab.|l]), where the uncertainty comes from 
our ignorance of the EOS and the mass of the star. 
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isGRID 


FIG. 25: Comparing the spin computation according to 
Eq. ( |C3| | (jfit) with the spin output by SGRID (j'sgrid) for 
SLy (red), ALF2 (blue), MS lb (green), and r2 (black). Ab¬ 
solute values (upper panel) and fractional residuals (bottom 
panel). 


Considering such systems, it is interesting to estimate 
the maximum spin that can be achieved with SGRID. 
For this purpose we consider once more the SLy, ALF2, 


and MSlb EOSs and single star configurations. The it¬ 
eration process to achieve high spins for realistic EOSs 
is as follows. We start with w = 0 (zero spin, i.e., a 
solution of the TOV-equation) and increase the angular 
velocity in steps of Aw = 0.005 up to w = 0.01 and in 
smaller steps of Aw = 0.0025 up to w = 0.0275. Note 
that higher spins could be computed for lower resolu¬ 
tions, but for higher resolutions the iteration procedure 
fails. Depending on the EOS, maximum dimensionless 
spins of 0.5 < jmax < 0.7 can be obtained (shown in the 


upper panel of Eig. 26). The middle panel shows the mass 
shedding parameter^ 


M 


X = 


^r^leq 
I pole 


(C4) 


which measures the deformation of the neutron star 
(caused by its rotation) according to the derivative of 
the enthalpy at the star’s surface parallel or perpendic¬ 
ular to the symmetry axis. Note that in binaries i9r^|eq 
is evaluated along the line connecting the two stars’ cen¬ 
ters. The lower panel shows the L^-norm Hamiltonian 
constraint in the domain covering the neutron star up to 
^max = 0.35, i.e., including the star’s surface, which is 
the most problematic region. The Hamiltonian constraint 
grows for angular velocities w > 0.02. We suggest that 
this is related to the deformed shape of the neutron star, 
indicated by the decreasing y. 
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